The overall flowchart of this study was shown in Figure 1 . Firstly, we screened out the differentially expressed genes (DEGs) between tumor and normal tissues in the TCGA-COAD cohort, followed by intersecting with LLPS-related genes obtained from DrLLPS database. Then, prognostic LLPS-related DEGs were identified by univariate Cox regression. Based on the expression profiles of these genes, unsupervised clustering analysis was performed to identify different LLPS patterns in colon cancer patients. Subsequently, we investigated the heterogeneity of two LLPS subtypes. In addition, the Least Absolute Shrinkage Selection Operator (LASSO) Cox regression is used to construct LLPS-related gene signature to calculate the LRRS. The robustness of signature is evaluated by multiple dimensions.

RNA sequencing gene expression profile (including count and TPM value), somatic mutation and clinical information of TCGA-COAD cohort were downloaded from the Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). The count value was used to identify the DEGs between tumor and normal tissues via “limma” R package (18). The normalized series matrix file of GSE39582 was directly downloaded from GEO database (https://www.ncbi.nlm.nih.gov/geo/) (19). The “ComBat” function of the “sva” R package was used to correct the batch effects of non-biological technical biases of TCGA-COAD TPM values and GSE39582 expression data. Patients without survival data were excluded from further analyses (20). Following the exclusion of normal tissue samples and data lacking overall survival (OS) information, the subsequent analysis was conducted on 448 samples from the TCGA-COAD cohort and 562 samples from the GSE39582 cohort. Data were analyzed with R (version 4.1.0).

A total of 4494 LLPS-related genes, of which 90 were scaffolds (2.00%), 487 were regulators (10.84%), and 3917 were clients (87.16%) ( Supplementary Table S1 ) in Homo sapiens involving 36 condensates were obtained from the DrLLPS database (http://llps.biocuckoo.cn/), which is a comprehensive database containing 437887 LLPS related proteins from 164 eukaryotes (21).

Firstly, the empirical Bayesian method of “limma” R package was used to identify DEGs between tumor and normal tissues in the TCGA-COAD cohort based on the screening criteria of P<0.05 and |log2FC|>1 (22). Then, intersect with 4494 LLPS-related genes, followed by univariate Cox regression analysis. A total of 253 prognostic LLPS-related DEGs were identified. Based on the expression of 253 prognostic LLPS-related DEGs mentioned above, unsupervised clustering analysis was performed using the “ConsensuClusterPlus” package to identify different subtypes of LLPS in colon cancer patients (23).

GSVA is a non-parametric, unsupervised method for estimating variation of gene set enrichment through the samples of an expression data set (24). To investigate the differences in Kyoto Encyclopedia of Genes and Genomes (KEGG) signaling pathways among different LLPS subtypes, we conducted GSVA enrichment analysis using the “GSVA” R package (24). The gene set “c2. cp. kegg. v2023.1. Hs. symbols” was download from the MSigDB database (https://www.gsea-msigdb.org/gsea/msigdb). Adjusted P-value less than 0.05 is considered statistically significant.

Gene Ontology (GO) annotates genes to three categories including biological processes, molecular functions, and cellular components (25). GO enrichment analysis was performed using the “enrichGO” function of the “clusterProfiler” R package (26).

Single Sample Gene Set Enrichment Analysis (ssGSEA), which is an extension of Gene Set Enrichment Analysis (GSEA), is used to assess the relative abundance of various immune cell infiltrations and immune functions in patients with colon cancer. Two gene sets, containing 23 and 29 types of immune cell markers respectively ( Supplementary Tables S2 , S3 ), are employed to evaluate the tumor microenvironment and immune functions of different LLPS subtypes in colon cancer (27, 28). The immune scores, stromal scores, ESTIMATE scores, and tumor purity of colon cancer patients was calculated by the “estimate” R package (29). The xCELL, TIMER, quanTIseq, EPIC, ConsensusTME, and ABIS methods from the “immunedeconv” R package are utilized to evaluate the correlation between the LRRS and immune cell infiltration (30). Expression levels of 45 immune checkpoint-related genes were analyzed between LLPS clusters (31).

The Tumor Immune Dysfunction and Exclusion (TIDE) algorithm (http://tide.dfci.harvard.edu/) was used to evaluate the tumor immune scape potential of colon cancer patients from their expression profiles (32, 33). Patients with lower TIDE scores was more likely to show stronger responses to immune therapy. The immunophenotypic score (IPS) data was download from The Cancer Immunome Atlas (TCIA) database (https://tcia.at/home). Then compared the IPS between different LLPS clusters to evaluate the responsiveness to anti-PD-1 or anti-CTLA-4 therapy (34).

To construct a LLPS-related gene signature, 253 prognostic LLPS-related DEGs were used to build a LASSO Cox regression model. LRRS was calculated as follows:

where Coef_i and Exp_i represent the LASSO coefficient and the corresponding gene expression level, respectively. Patients with colon cancer were stratified into low and high LRRS groups according to the median value of LRRS. Univariate and multivariate Cox regression analyses were conducted on LRRS in conjunction with clinical characteristics to identify independent prognostic factors. A nomogram was constructed utilizing these independent prognostic factors. Receiver operating characteristic (ROC) curve analysis were applied to assess the accuracy of nomogram, LRRS, and clinical characteristics in predicting OS of colon cancer patients. The concordance index (C-index) was employed to assess the prognostic capability of nomogram, LRRS, and clinical characteristics. The calibration curves were used to evaluate the precision of the nomogram in terms of the agreement between the observed and predicted OS outcomes at the 1st, 3rd, and 5th years.

The mutation profiles and frequencies were visualized using the “maftools” R package. The tumor mutation burden (TMB) was defined as the number of mutations per megabase (mut/Mb) (35). The copy number variation (CNV) data was acquired from UCSC Xena (https://xena.ucsc.edu/) database. Then identify the copy number amplifications or deletions of model gene in the cohort.

HCT-8 cells were purchased from Wuhan Pricella Biotechnology Co., Ltd. The cells were cultured in RPMI-1640 medium (Gibco, USA) supplement with 10% Fetal Bovine Serum (FBS, Invitrogen Corporation, USA) and 1% penicillin/streptomycin (Gibco, USA) and incubated in humidified incubator containing 5% CO2 at 37 °C.

Cells (5 × 10⁵ cells/well) were seeded in the six-well plate. After cell attachment, lentivirus carrying overexpression plasmid was added to the culture medium, supplemented with polybrene to reach a final concentration of 1 µg/mL. The lentivirus was removed and replaced with normal growth medium at 12 h transfection. 48 h after transfection, 1 μg/mL puromycin was added to screen for stable cells for one week.

Total RNA was extracted using RNA Extraction Reagent TRIzol (Invitrogen, USA). Then RNA concentration was determined with a spectrophotometer. The RNA was further reverse-transcribed to cDNA using cDNA synthesis kit (Bio-Rad, USA). cDNA was amplified with SYBR Green Supermix (Bio-Rad, USA). The 2^−ΔΔCt method (ΔΔCT = ΔCt control−ΔCt sample) was used to calculate the amplification fold. Primer sequences are shown in Table 1 .

The cells were seeded at 1×10⁵ cells/100 μl in 96 well plates (100 μL/well) and incubated for 24 hours. Add 10 μL of CCK-8 solution (Dojindo, China) to each well and incubate at 37°C for 1-4 hours. The absorbance was determined at 450 nm using the absorbance microplate reader (Bio-Rad, USA).

Cell suspension was diluted to the final concentration of 10³ cells/mL. Then the cells were seeded into six-well plates (1 mL/well) and incubate for 10 days. After the clones were formed, the cells were fixed with 4% paraformaldehyde for 30 min and stained with crystal violet staining solution for 20 min. Finally, the cells were washed with PBS and photographed. The experiment was repeated at least three times.

5x10⁵ cells were seeded in each well of a six-well plate and when the cells were grown to 90% confluency, a straight line was scratched across the cell monolayer with a 200 μL pipette tip. Then, add medium containing 1% serum to each well after washed with PBS. Cell migration was observed and photographed at 0 h and 24 h under the microscope. The experiment was repeated at least three times.

2 × 10⁴ cells were added to the upper chamber with serum-free medium (Gibco, USA). 600 μL of complete medium was added to the lower chamber. After 24 hours of incubation, the cells in the upper chamber were removed with cotton swabs, and the cells on the lower surface of the chamber were fixed with 4% paraformaldehyde for 30 minutes, and then stained with crystal violet for 30 minutes. Finally, five visual fields were randomly selected to be photographed with the microscope (Olympus, Japan).

The Kruskal-Wilcoxon test was used to compare statistical differences between two groups. While statistical differences among more than two groups was compared using the Kruskal-Wallis test. The correlation between expression of model genes was assessed through the Spearman correlation analysis. The survival curves for the prognostic analysis were generated via the Kaplan-Meier method and the significance of differences were identified by log-rank tests. All statistical p value were two-side, with p < 0.05 as statistically significance. All data processing was done in R 4.1.0 software.

To identify the LLPS-related DEGs, we conducted the following screening. Firstly, according to the screening criteria of adj. p<0.05 and | log2FC |>1, 8002 DEGs between tumor samples and normal samples were obtained from the TCGA-COAD cohort. Among them, 3584 were upregulated and 4418 were downregulated in tumor samples compared with normal samples ( Supplementary Table S4 ). Subsequently, after intersecting with 4494 LLPS-related genes, 980 LLPS-related DEGs were identified, of which 444 were upregulated and 536 were downregulated in tumor samples ( Supplementary Table S5 ). Then, these LLPS-related DEGs were subjected to univariate Cox regression analysis to identify 253 prognostic LLPS-related DEGs ( Figure 2A ), of which 3 were scaffolds, 26 were regulators, and 224 were clients ( Supplementary Table S6 ).

Based on the expression levels of the 253 prognostic LLPS-related DEGs mentioned above, we performed unsupervised clustering analysis on the merged TCGA-COAD and GSE39582 cohorts, and ultimately divided 1010 colon cancer patients into two subtypes: LLPS cluster A (n=744) and cluster B (n=266) ( Figure 2B ). The significant differences in the expression of 253 prognosis LLPS-related DEGs between the two subtypes were observed in the heatmap ( Figure 2C ). Prognostic analysis indicated a significant survival difference between the two subtypes of LLPS. Cluster A had a better overall survival outcome than cluster B ( Figure 2D ).

To explore the potential molecular mechanisms underlying the LLPS subtype of colon cancer, we performed GSVA to evaluate the differential KEGG gene sets between the two subtypes. The results showed that cluster A exhibited associations with cell cycle, DNA replication and mismatch repair, whereas cluster B demonstrated associations with the MAPK signaling pathway and JAK-STAT signaling pathway ( Figure 2E ). Additionally, we performed GO enrichment analyses on 253 prognosis LLPS-related DEGs ( Figure 2F ).

We attempted to compare the clinicopathological characteristics of colon cancer patients in two different LLPS clusters. Compared to cluster B, cluster A had more patients with age greater than 65 years and higher proportion of stage I and stage II, which may explain why the overall survival outcome of cluster A is better. No statistically significant variances were observed between the two clusters in terms of gender and MSI status distribution ( Figure 3 ).

Subsequently, we conducted a comparative analysis of immune cell infiltration among the two subtypes of LLPS. Results derived from two distinct gene sets consistently indicated that cluster B exhibited a higher degree of immune cell infiltration ( Figures 4A, B ). Employing the ESTIMATE algorithm, we predicted the abundance of stromal and immune cells across different subtypes. The analysis revealed that cluster A exhibited lower ESTIMATE, immune, and stromal scores compared to cluster B, indirectly reflecting a higher tumor purity in cluster A ( Figures 4C, D ). Then, we compared the immune functions between the two LLPS subtypes. Notably, cluster B demonstrated significantly elevated functions, particularly in aspects such as immune check-point ( Figure 4E ). A detailed differential analysis highlighted that the expression levels of key immune checkpoint molecules, including CD274 (PD-L1), CTLA4, LAG3, and TIGIT, were markedly higher in cluster B compared to cluster A ( Figure 4F ). Finally, we compared the IPS between the two LLPS clusters. Cluster A demonstrated a significantly elevated IPS, suggesting a more promising immunotherapeutic response potential ( Figure 4G ). TIDE analysis indicated that cluster A had a lower TIDE score, further corroborating its enhanced likelihood of responding favorably to immunotherapeutic interventions ( Figure 4H ).

Next, 253 prognostic LLPS-related DEGs were subjected to LASSO Cox regression to construct a LLPS-related gene signature for predicting the prognosis of colon cancer. According to the ratio of 1:1, 1010 colon cancer patients were randomly divided into a train group and a test group. The detailed clinical information of the train group, test group and total group is shown in Table 2 . Ultimately, 14 genes were identified ( Figures 5A, B ), comprising 7 protective genes and 7 risk genes, with their respective risk coefficients detailed in Table 3 . The LRRS of each colon cancer patient can be calculated by summing the product of the expression levels of each gene and the corresponding risk coefficients. The area under the curve (AUC) values were evaluated by ROC curve. The LRRS had the highest AUC value in the 3rd year, reaching 0.730. Additionally, AUC values of 0.677 and 0.697 were observed in the first and fifth years, respectively ( Figure 5C ). According to the median value of LRRS, colon cancer patients in the train and test group were divided into high and low risk groups, respectively. The distribution of LRRS, survival time, status, and expression heatmaps of risk model genes among high- and low-risk patients in the total, train, and test groups were displayed ( Figures 5D–F ). The prognostic analysis of the overall survival time of CRC patients revealed that those with high LLPS scores exhibited poorer prognoses in all the three groups ( Figure 5G ). Furthermore, a comparative analysis of progression free survival in the TCGA-COAD cohort ( Figure 5H ) and recurrence-free survival in the GSE39582 cohort ( Figure 5I ) demonstrated that individuals with higher LRRS experienced inferior prognoses in the total, train and test groups.

We want to build a survival prediction model for colon cancer patients that can be applied in clinical practice. To this end, we first conducted univariate and multivariate Cox regression analyses of LRRS and clinical information with overall survival. The results indicated that the LRRS is an independent prognostic factor for OS. In univariate Cox regression analyses, the hazard ratio (HR) of LRRS was 1.154 with a 95% confidence interval (CI) of 1.122-1.187 (p<0.001, Figure 6A ). In multivariate Cox regression analyses, the HR of LRRS was 1.127 with a 95% CI of 1.087-1.169 (p<0.001, Figure 6B ). In addition, age and TNM stage were also independent prognostic factors. Subsequently, a nomogram was developed by integrating age, TNM stage, and LLPS risk ( Figure 6C ). By calculating the score of each variable, the 1-, 3-, and 5-years survival of colon cancer patients can easily estimate by drawing a vertical line. The ROC curve demonstrated that the LRRS had excellent accuracy in predicting OS, with an AUC of 0.680, surpassing the predictive ability of any other clinical feature. Furthermore, the nomogram showed an even higher AUC of 0.792 ( Figure 6D ). Meanwhile, the concordance index indicated that the nomogram had the highest predictive accuracy, followed by LRRS ( Figure 6E ). The calibration curves also illustrate that there is a good agreement between the observed and predicted OS outcomes at the 1st, 3rd, and 5th years for the nomogram ( Figure 6F ).

We conducted a comparison of the LRRS among various LLPS clusters, revealing that cluster B exhibited a significantly higher LRRS compared to cluster A ( Figure 7A ). The relationship between LLPS clusters, LRRS, clinical stage, and survival status is shown in Figure 7B . Notably, a majority of patients within Cluster B were categorized into the high LRRS group, consistent with the results in Figure 7A . Then, the Wilcoxon test was used to calculate the correlation between LRRS and clinicopathological features including age, gender, clinical stage, and survival status. The results showed that high LRRS was associated with high clinical stage and poor prognosis, but not with age and gender. The LRRS gradually increased from stage I to stage IV ( Figure 7C ).

Subsequently, employing a panel of methodologies including xCELL, TIMER, quanTIseq, EPIC, ConsensusTME, and ABIS, we investigated the correlation between LRRS and immune cell infiltration. The findings indicated that T cell NK, macrophage, dendritic cell, cancer associated fibroblast, and monocyte were positively correlated with LRRS, while basophil and T cell CD4+ memory were negatively correlated with LRRS ( Figure 7D ). An immune functional correlation analysis revealed that LRRS was positively correlated with immune functions such as immune checkpoint and T cell co-inhibition ( Figure 7E ). Consistently, LRRS was positively correlated with the expression of immune inhibitory molecules such as CD274, LAIR1, and NRP1 ( Figure 7F ). In our final analysis, we assessed the relationship between LRRS and comprehensive immunological scores, including immune, stromal, ESTIMATE scores, and tumor purity. The results demonstrated that LRRS was positively correlated with immune, stromal, and ESTIMATE scores, but negatively correlated with tumor purity ( Figure 7G ). In summary, these analyses underscore that while the high LRRS group shows enhanced immune cell infiltration, it paradoxically exhibits a higher state of immune suppression.

Next, we analyzed the top 20 mutated genes in the two different groups ( Figure 8A ). The gene mutation frequencies were found to be similar between the two groups. APC, TP53, and TTN are the three genes with the highest mutation frequency. Additionally, the TMB was evaluated between the high and low LRRS groups, revealing no significant difference in TMB levels ( Figure 8B ). However, survival analysis showed that the prognosis of the high TMB group was significantly worse than that of the low TMB group ( Figure 8C ). Therefore, we wonder whether it is possible to combine TMB and LRRS to stratify patients, and the results showed that this combination could better predict patient prognosis. The group with high TMB and high LRRS had the worst prognosis, followed by the group with low TMB and high LRRS ( Figure 8D ).

A total of 14 genes were utilized in the construction of the LLPS-related gene signature. Among these genes, AGAP3, CDK2, DMKN, PRMT1, PSMA7, POU4F1, RAB15 and SNAI1 were up-regulated in tumor tissue compared with normal tissue in TCGA-COAD cohort, while the remaining genes were down-regulated in tumor tissue ( Figure 9A ). Hazard ratios (95% CI) of these genes were calculated by univariate Cox hazard analysis ( Figure 9B ). In terms of gene types, these 14 genes consisted of 11 clients, 2 regulators, and 1 scaffold. Except for the scaffold SYN2, there existed an expression correlation among other genes ( Figure 9C ). Next, we summarized the incidence of somatic mutations of 14 LLPS-related genes in TCGA-COAD cohort. Among 399 tumor samples, only 64 (16.04%) exhibited mutations in the model genes. The gene with the highest mutation frequency was NRG1, but only 5%, indicating that these LLPS-related model genes are conserved in the progression of colon cancer ( Figure 9D ). The statistical results of CNV alteration frequency showed that the CNV alteration was prevalent in 14 model genes. Specifically, copy number amplification was more pronounced in CACNB1, CDK2, DMKN, and SNAI1, while CLMN, NRG1, OGDHL, PRMT1, and SYN2 mainly exhibit copy number deletion, consistent with their mRNA expression ( Figure 9E ). The location of CNV changes of these genes on chromosomes was shown in Figure 9F .

Given that POU4F1 is a regulator which may play a key role in the LLPS process, we conducted an overexpression study of POU4F1 in HCT-8 cells to elucidate its effects on cell proliferation and migration. POU4F1 overexpression were achieved by lentivirus transfection ( Figure 10A ). Further investigations were conducted to determine the proliferation and migration of HCT-8 cells. The CCK8 assay showed cell proliferation increased after POU4F1 overexpression ( Figure 10B ). Meanwhile, an elevated number of clones were observed followed by POU4F1 upregulation ( Figure 10C ). Subsequently, wound healing assay was performed and the results showed that the overexpression of POU4F1 resulted in a significant enhancement of the cell migration ( Figure 10D ). Transwell migration assay also showed that the number of migrated cells in the POU4F1-overexpressing group was greater than that in the control group ( Figure 10E ). The findings above suggested that POU4F1 facilitated both proliferation and migration of HCT-8 cells.