The tumor ploidyscore of each tumor were obtained from the previous study [19]. Pearson correlation analysis between BEX4 expression and ploidy score was performed to identify significant correlations. The GEPIA (https://genomics.senescence.info/cells/) and HUMANBASE (https://hb.flatironinstitute.org/) databases provided a list of 178 BEX4-related genes. Aneuploidy scores were collected from the cBioPortal database (https://www.cbioportal.org/). Machine learning models (random-forest (RF), support vector machine (SVM)) were applied to analyze gene matrices and identify aneuploid characteristic genes (ACGs), with model fitting assessed via residual analysis and receiver operating characteristic (ROC) curve visualization using the Caret (version 6.0), ggplot2 (version 3.5.1), randomForest (version 3.3.1), kernlab (version 0.9), and pROC (version 1.18.5) R packages.

As a classical ensemble algorithm, Random-Forest (RF) has superior regression and classification capabilities by constructing multiple decision trees and forming forests. And We averaged all aneuploidyscore of samples from different cancers and grouped them into high and low groups. We inputted a matrix containing 178 gene expressions in 9 cancers. The samples with low aneuploidyscore were set as the control group, and the samples with high aneuploidyscore were set as the training group. With support from ConsensusPathDB (http://cpdb.molgen.mpg.de/), we were able to perform enrichment analysis of aneuploidy characteristic Hub-genes (AHGs) at the gene functional level and study the cell components that are mainly affected by its expression.

We extracted 30 ACGs in differential tumors and plotted Nomogram. At the same time, as an important part of the prediction model, we used the Hosmer-Lemeshow test to test the gap between the actual value and the predicted value and made the continuous calibration curve. The decision curve analysis (DCA) and clinical impact curve (CIC) were used to further evaluate the clinical predictive effect of the model. In this process, we use the R package: ‘rms’ (version 1.1), ‘rmda’ (version 1.6).

Using 30 ACGs as input variables for kidney chromophobe (KICH), thyroid carcinoma (THCA), and thymoma (THYM), consistency analysis was performed using the ConsensusClusterPlus R package(version 1.66.0), with cluster numbers set to 2 and 100 iterations extracting 80% of the sample. Visualization was done via pheatmap (version 1.0.12) and ggplot2.

The Mfuzz package clustered Mfuzz expression patterns for KICH, THCA, and THYM samples from the TCGA database. ssGSEA scores for different clustering modules and tumor patients were calculated, identifying the gene module most correlated with BEX4.

We obtained four important gene sets from Miranda et al.’s study [21] that can represent the evaluation of tumor stemness and calculated the corresponding ssGSEA-stemness-scores for all samples based on the expression matrix, including: Cell cycle score, Cell Transcription Score, Metabolism Synthesis Score and Signal recognition and transduction score. Next, we collected aneuploidy scores from the CI database, and comprehensively constructed a sample matrix containing 55 scores (including 50 cluster pattern expression scores, 4 stemness scores and aneuploidy scores) for all samples in pan-cancer. And Pearson correlation analysis was performed between the two.

Using 30 ACGs as input variables, we loaded the R package ‘ConsensusClusterPlus’ (version 1.66.0) for consistency analysis, set the number of clusters to 2, repeat 100 times to extract 80% of the total sample in KICH, THCA and THYM. The clustering heat map is visualized by ‘pheatmap’. Gene expression heat map retained genes with variance above 0.1. The box line diagram is drawn by ‘ggplot2’. Then, we used non-parametric tests between different types to understand the differential expression of AHGs.

The aneuploidy scores of all collected TCGA samples were analyzed in the form of non-parametric tests between the two subtypes of different cancers to verify whether the sample typing using ACGs can distinguish the high aneuploidy score group and the low aneuploidy score group at the gene expression level. At the same time, the Principal Component Analysis (PCA) score was constructed as aneu-score in KICH, THCA and THYM using the characteristic genes screened by RF, and its correlation with the vest was analyzed.

Expression values of eight immune checkpoint-related genes (SIGLEC15, TIGIT, CD274, HAVCR2, PDCD1, CTLA4, LAG3, PDCD1LG2) were extracted and visualized using ggplot2 and pheatmap. The TIDE algorithm predicted Immune Checkpoint Blockade (ICB) responses, with immune scores calculated using MCP-counter, CIBERSORT, TIMER, and quanTIseq.

Several studies have demonstrated that the occurrence of polyploidy is a significant contributing factor to the development of Cell Stemness Score (CSCs), alongside the transformation of tissue stem cell [22, 23]. To determine the mRNAsi of the samples, the logistic regression machine learning algorithm (OCLR) [24, 25] was employed, which utilize the pan-cancer expression and clinical information files available on TCGA. When we obtained the mRNAsi scores of all samples, a difference analysis was performed between different subtypes to understand the association between the degree of aneuploidy acquisition and the degree of stemness acquisition. We employed Spearman correlation analysis mapping the dryness index onto the (0–1) interval.

HPA (https://www.proteinatlas.org/) is a human proteome atlas database that contains information on the distribution of proteins in human tissues and cells. With the help of HPA database, we collected immunohistochemical images of 16 tumor tissues. Meanwhile, HPA database also showed that it had significant prognostic efficacy only in pancreatic cancer and renal cancer (p 0.001). This study included 116 tumor samples (29 PAADs, 29 KICHs, 29 KIRCs, and 29 KIRPs) collected from the First Affiliated Hospital of Anhui Medical University from 2018 to 2023. Sections were taken from the TMA paraffin blocks, deparaffinized in xylene and rehydrated in a descending series of alcohol. After antigen retrieval by using citrate buffer, the sections were incubated in 3% H₂O₂ to block endogenous peroxidase activity. Then the sections were incubated with primary antibody (BEX4, Bioss, Beijing, China, bs-7978R, dilute 1:100) for 1 h at room temperature. After incubated with a horseradish peroxidase (HRP)-labeled polymer and 3,3-diaminobenzidine chromogen solution (Dako EnVision), the sections were counterstained with Harris hematoxylin (Solarbio® Life Sciences, Beijing, China). The immunohis-tochemical analysis of BEX4 in the sections was performed using a semiquantitative scoring system to distinguish the staining groups as the strong and weak groups. Cytoplasmic staining was positive, and immunopositive degree was divided into weak (1) and strong (2) according to staining intensity.

Western blotting was conducted as previously described [26]. cells were lysed in radio immunoprecipitation assay lysis (RIPA) buffer (P0013), and protein concentrations were measured using the bicinchoninic acid (BCA) assay. Proteins were separated via 10% or 15% sodium dodecyl sulfate–polyacrylamide gel electrophoresis (SDS-PAGE) and transferred to nitrocellulose membranes (Mil-lipore, Burlington, MA, USA). Following an overnight incubation with the primary antibody (BEX4, Bioss, bs-7978R, dilute 1:100) at 4 °C, the membranes were exposed to an HRP-conjugated secondary antibody (BEX4, Bioss, Cat No. SA00001-2). Protein complexes were visualized using enhanced chemiluminescence reagents (Millipore, USA).

Human renal cell carcinoma cell line A498 and pancreatic cancer cell line PANC-1 were purchased from Biyuntian Biotechnology Company (Shanghai, China). A498 or PANC-1 cells (2.5 $\times$ 10⁵/well) were seeded in the upper chamber of a 24-well transwell (Corning, NewYork, NY, USA) with DMEM lacking FBS. All cell lines were verified by short tandem repeat (STR) profiling and no mycoplasma contamination was detected. Cells were cultured at 37 °C with 5% CO₂. The lower chamber contained DMEM with 10% FBS. After 24 hours, migrated cells on the underside of the membrane were fixed in methanol, stained with 0.1% crystal violet, and imaged under a microscope (Leica, Wetzlar, Germany).

Cells (2 $\times$ 10³) were pre-seeded into 96-well plates. After incubation for 24, 48, and 72 hours, the medium was replaced with cell counting kit-8 (CCK-8) reagent (MCE, Shanghai, China), 100 µL per well. OD values at 450 nm were measured after 1 hour of incubation at 37 °C.

Statistical analyses were performed using R (version 4.3.1) (https://www.r-project.org/) and SPSS 27 (IBM Corp., Chicago, IL, USA). All relative risks were presented with a 95% confidence interval (p 0.05). Data were expressed as mean $\pm$ standard deviation (SD). Group differences were evaluated using Student’s T-test, with p 0.05 considered significant. Tumor abbreviations and their corresponding full names are listed in Supplementary Table 1.

After conducting correlation analysis with BEX4, a total of 178 genes were identified as being related. Among these genes, there were 9 significant tumors, namely THYM, KICH, OV, TGCT, PAAD, THCA, SARC, LIHC, and UCEC. These tumors were selected based on their strong correlation between ploidyscore and the expression of BEX4, as indicated in Supplementary Table 2. Further analysis using residual tests demonstrated that RF efficacy was more significant in all cancer types except UCEC, which was confirmed by the ROC curve (Fig. 3A–I and Supplementary Fig. 3A–I). Subsequently, RF or SVM was utilized to identify the respective aneuploid characteristic genes in the 9 cancer species (Supplementary Fig. 4J–O). Subsequently, we identified the recurrent presence of three genes, namely BEX4, CAMSAP2, and MARCKS, as AHGs through the intersection of aneuploidy characteristic genes across different tumors. To further investigate their significance, we selected KICH, THCA, and THYM for subsequent studies (Fig. 3J–L, respectively), considering multiple factors. Additionally, we developed nomogram prediction score utilizing the aneuploidy characteristic genes (Supplementary Fig. 4P–R). The effectiveness of the model was further confirmed through decision curve analysis (DCA) and clinical impact curve (CIC) in KICH (Fig. 3M), THCA (Fig. 3N) and THYM (Fig. 3O).

Fig. 3.

The efficiency test and application of machine learning. (A–I) The curves shows the performance test of random forest and support vector machine using residual analysis in 9 cancer species. (J–L) The recurrent presence of three genes, BEX4, CAMSAP2, and MARCKS, identified as Aneuploidy-Related Characteristic Genes (AHGs) through random forest (RF) or support vector machine (SVM) analysis, is illustrated in KICH, THCA, and THYM. These tumors were selected for further studies due to their strong correlation with BEX4 expression and aneuploidy characteristics. (M–O) DCA (decision curve analysis) and CIC (clinical impact curve) were used to further test the effect of feature gene screening in three kinds of tumors. KICH, kidney chromophobe cell carcinoma; THCA, thyroid carcinoma; THYM, thymoma.

The enrichment results, as facilitated by ConsensusPathDB, encompassed various cellular components and molecular functions related to the microtubule cytoskeleton, including cytoskeletal part, centrosome, cytoskeleton, microtubule organizing center, spindle pole, cytoskeletal protein binding, calmodulin binding, tubulin binding, spindle, microtubule, intracellular non-membrane-bounded organelle, supramolecular fiber organization, and polymeric cytoskeletal fiber (Supplementary Table 3). The gene function enrichment analysis revealed that three AHGs had an impact on the biological process of cytoskeleton protein binding and the cytoskeleton in cell components. Additionally, BEX4 and CAMSPA2 were found to influence the synthesis of polymer cytoskeleton fibers. Furthermore, the analysis of cell components indicated that the three AHGs could affect the microtubule cytoskeleton, with MARCKS and CAMSAP2 affecting both the microtubule organization center and centrosome. Moreover, BEX4 and CAMSAP2 were found to have a combined effect on microtubules. Lastly, BEX4 and CAMSPA2 were found to influence the molecular function of microtubule protein binding.

We found 4 miRNAs strongly correlated with BEX4, CAMSAP2, and MARCKS in the ENCORI database: hsa-miR-425-5p, hsa-miR-200c-3p, hsa-miR-200b-3p, and has-miR-425. From these miRNAs, we identified 6 related genes and 6 lncRNAs. Using these 3 genes, we found more miRNAs and created the LncRNA-miRNA-Gene Co-expression network (Fig. 4A,B). After identifying ACGs, we established the ceRNA and gene network, which demonstrated the mutual regulation of AHGs in tumor tissues (Fig. 4C). After excluding miRNAs with less than 30% expressed individuals, TCGA-based survival analysis identified significant risk factors in different cancers: hsa-miR-200b-3p in THYM, hsa-miR-200c-3p in LIHC and THYM, hsa-miR-425-5p in KICH (but protective in OV), and hsa-miR-429in THYM. Then, to explore the prognosis value of those miRNAs and lncRNAs, we employed survival analysis and KM curve to investigate the survival status of samples with different expression of miRNAs or lncRNAs respectively (Fig. 4D–I).

Fig. 4.

Construction of ceRNA-genes Co-expression network and prognostic significance of miRNA in KICH, THCA and THYM. (A) The comprehensive regulatory network suggested miRNAs that may regulate the expression of ACGs. The yellow triangle module represented lncRNA, the green rectangle module represented miRNA, and the green polygon module represents gene. The network demonstrates the intricate interactions and co-expression patterns among long non-coding RNAs (lncRNAs), microRNAs (miRNAs), and genes, highlighting their potential regulatory roles in cancer biology. (B) Through the extraction of the integrated network, we obtained a close and multi-layer accurate ceRNA-genes network (C) The ceRNA-genes network constructed by AHGs. (D–I) Survival analysis identified significant risk factors in various cancers, with miR-200b-3p in THYM, miR-200c-3p in LIHC and THYM, miR-425-5p in KICH and OV, and miR-429 in THYM, respectively. ACGs, aneuploid characteristic genes; AHGs, aneuploidy-related characteristic genes. LIHC, liver hepatocellular carcinoma; OV, ovarian serous cystadenocarcinoma.

ACGs identified two molecular subtypes in KICH, THCA, and THYM (Fig. 5A-1,2; B-1,2; C-1,2). The aneuploid score from the Cbioportal database and a difference test revealed significant differences in aneuploidy between the subtypes (Fig. 5A-3, B-3, C-3). Notably, in KICH and THCA, the subtype with a higher aneuploidy score had a lower tumor cell stemness. At the same time, the difference between the two subtypes of AHGs in the three cancer species was not consistent (Fig. 5A-6,B-6,C-6). The correlation analysis results showed a correlation coefficient greater than 0.3 (Fig. 5A-7,B-7) between the Aneu-score based on characteristic genes and aneuploidies. In THYM, the correlation coefficient reached 0.6 (Fig. 5C-7). Additionally, in pan-cancer, the immune checkpoint response indicated differential responses between the two subtypes in their respective tumors (Fig. 5A-8,B-8,C-8). In THYM, the TIDE SCORE difference indicated a significant difference in the effectiveness of immune checkpoint blockade therapy (ICB) between the two patient types, with C1 being superior to C2 (Fig. 5C-9). Additionally, CAMSAP2 expression was significantly higher in the immunotherapy response group compared to the non-response group (Supplementary Fig. 8G). The analysis of immune infiltration scores revealed that Macrophage and Neutrophil were the main differing immune cells (Fig. 5A-10,B-10,C-10, and Supplementary Fig. 8A–C).

Fig. 5.

Identification of Molecular Subtypes and Analysis of Aneuploidy, Tumor Cell Stemness, and Immune Response in KICH, THCA, and THYM. (A–C) The three modules represented Two molecular subtypes were identified in KICH, THCA, and THYM using Aneuploidy-Related Characteristic Genes (ACGs). (A-1,2; B-1,2; C-1,2) All the heatmaps showed the results of cluster in KICH, THCA and THYM. (A-3,B-3,C-3) MRNAsi scores were collected and analyzed between the two subtypes of different tumors. (A-4,B-4,C-4) The difference violin plot showed the difference in aneuploidy scores between different subtypes. (A-5,B-5,C-5) The Kaplan-Meie (KM) curve compared and visualized the survival information of patients with different subtypes in three cancer types. (A-6,B-6,C-6) Multiple sets of box plots in three cancer species visualized the differential expression of AHGs between subtypes. (A-7,B-7,C-7) Correlation analysis between the Aneu-score, based on characteristic genes, and aneuploidies demonstrated a correlation coefficient greater than 0.3 in KICH and THCA, while in THYM, the correlation coefficient reached 0.6, indicating a strong association between aneuploidy and the expression of ACGs. (A-8,B-8,C-8) The expression of immune checkpoint in different subtypes was significantly different and there was a great difference between different tumors. (A-9,B-9,C-9) checkpointblockadetherapy (ICB) treatment varied among different tumor subtypes. The above pictures were the statistical data of immune response of samples in different subtypes, and the following pictures were the box diagram of the difference results of treatment response scores. (A-10,B-10,C-10) The immune infiltration score analysis highlighted macrophages and neutrophils as the main differing immune cells between the subtypes, suggesting their potential role in the immune response and tumor microenvironment in these cancers. *p 0.05, **p 0.01, ***p 0.001, ****p 0.0001.

In our study, we performed clustering based on the expression level of BEX4 who Mfuzz expression pattern, and successfully obtained 50 clustering results in KICH (Supplementary Fig. 5A), THCA (Supplementary Fig. 6A), and THYM (Supplementary Fig. 7A), respectively (Supplementary Table 4A: Gene lists in different cluster pattern.xls). Combined with the ssGSEA score and expression characteristics of the clustering module, we found that in pan-cancer, the score obtained by clustering according to the BEX4 expression pattern (Supplementary Table 4B: Correlation between the expression of BEX4 and ssGSEA score.xls) can effectively distinguish the degree of aneuploidy in the tissue (Supplementary Fig. 5C,6C,7C). At the same time, there are many expression pattern clustering gene sets that are significantly correlated with the expression trend of BEX4 in each cancer species (Supplementary Fig. 5B,6B,7B).

In KICH, there was no significant correlation between aneuploidy score and stemness score. However, in the four scores related to stemness, Cell Transcription Score was significantly positively correlated with Cell cycle score and Metabolism Synthesis Score, and Cell cycle score was also significantly positively correlated with metabolism. The expression score of clusters 24, 43 clustered gene set was significantly positively correlated with the score of stemness thought and was significantly correlated with the expression of BEX4. Through the enrichment analysis of clusters 24, 43, we observed that the genes in the cluster were involved in tumor metabolism and infection-related pathways. Then, for aneuploidy score, clusters 17, 19 was positively correlated with GSEA score, and showed a significant correlation with BEX4 expression, and was enriched in many metabolic related pathways (Supplementary Fig. 5D,E).

In THCA, there is a correlation between aneuploidy score and signal recognition pathway. In the stemness score, Cell Transcription Score was significantly positively correlated with Cell cycle score and Metabolism Synthesis Score, and Cell cycle score was also significantly positively correlated with metabolism. And Signal recognition and transduction score also showed a significant correlation with the other three scores. The gene expression scores in clusters 1, 13, 23 showed a significant positive correlation with the four scores of stemness, and a strong correlation with BEX4 expression. The genes in clusters 2, 6, 8, 18, 21, 26, 32, 50 showed a significant negative correlation with stemness score and a strong correlation with BEX4 expression. Cluster 11 showed a strong correlation with aneuploidy score and BEX4 expression. Enrichment results showed that cluster 11 was mainly enriched in related biological processes such as metabolism and virus, such as: Arginine and proline metabolism; legionellosis; hepatitis and so on (Supplementary Table 4C: Samples with different scores in clsuters.xls). There was a strong positive correlation between clusters 50 and 26, 32. The genes in the three clusters were mainly enriched in metabolic-related pathways and biosynthesis (Supplementary Fig. 6D,E).

In THYM, we observed a strong positive correlation between cell cycle and cell transcription, and a negative correlation between cell cycle and metabolism (Supplementary Fig. 7D,E).

In the process of cancer progression, there is a gradual loss of differentiated phenotypes and a gain of progenitor/stem cell-like characteristics. Through the calculation of the mRNAsi score, it was observed that there were predominantly negative correlations in most tumors. Notably, the correlation coefficient reached –0.75 in ESCA, –0.57 in THYM, –0.21 in THCA, and –0.68 in STAD. These strong correlation coefficients indicate a close association between BEX4 and the stem cell score across multiple tumors (Supplementary Fig. 3A–Y). Furthermore, in KICH, THCA, and THYM, calmodulin regulated spectrin associated protein 2 (CAMSAP2), and myristoylated alanine rich protein kinase C substrate (MARCKS) exhibited a significant negative correlation with the acquisition of dryness (Supplementary Fig. 3Z–AD).

In KICH, only LncRNA (BOLA3-AS1) showed significant high expression in low aneuploid subtype (Supplementary Fig. 9E). However, in THCA, miRNAs (miR-429 and miR-200b-3p) (Supplementary Fig. 9K,L) and LncRNAs (BOLA3-AS1, PCAT19 and RNF216P1) (Supplementary Fig. 9M–O) showed high expression in high aneuploid subtypes. Finally, in THYM, LncRNAs (BOLA3-AS1 and RNF216P1) (Supplementary Fig. 9U,W) were significantly down-regulated in high aneuploidy subtypes. Other miRNAs (miR-200c-3p, miR-425-5p and miR-429) (Supplementary Fig. 9Q–S) and LncRNA (PCAT19) (Supplementary Fig. 9V) were differentially expressed in high aneuploidy subtypes. And the differential expression trend of miRNAs (miR-425-5p, miR-429 and miR-200c-3p) between the two subtypes of THYM was opposite to the expression trend of AHGs.

To evaluate the protein expression of BEX4, we used the HPA database (https://www.proteinatlas.org/) to extract immunohistochemical image and found that the protein expression of BEX4 in tumor tissues of 16 cancer types was significantly lower than that in normal tissues (Supplementary Fig. 10A). Meanwhile, we also found that patients with high level of BEX4 had a good prognosis in pancreatic cancer and renal cancer in HPA database KM plotter database (Supplementary Fig. 10B,C).

To investigate the influence of BEX4 expression on the prognosis of patients in PAAD, KICH, KIRC, and KIRP, level of BEX4 was conducted by immunohistochemistry (IHC) and patients were divided into two groups with different level of BEX4 (Fig. 6A). In PAAD, KICH, KIRC, and KIRP, we found that patients with high level of BEX4 had great prognosis than patients with low level of BEX4, which indicated that high expression of BEX4 was a protective factor (Fig. 6B–E). These results were consistent with the information provided by the HPA database (Supplementary Fig. 10B,C). Then, to confirm the influence of BEX4 expression on the renal cancer and pancreatic cancer, we constructed BEX4 overexpression in A498 cells (A498 OE) and PANC-1 cells (PANC-1 OE), respectively (Fig. 6F,G). Subsequently, we found that overexpression of BEX4 alleviated motility and proliferation of A498 cells and PANC-1 cells (Fig. 6H–J), respectively.

Fig. 6.

The influence and prognostic value of BEX4 in pancreatic cancer and renal cell carcinoma. (A) Immunohistochemistry (IHC) showed the low and high expression of BEX4 in pancreatic adenocarcinoma (PAAD), KICH, kidney renal clear cell carcinoma (KIRC) and kidney renal papillary cell carcinoma (KIRP). Scale bar = 200 µm and 400 µm . (B–E) KM survival curves showed patients with high expression of BEX4 has great prognosis in PAAD (B), KICH (C), KIRC (D) and KIRP (E), respectively. (F,G) Western blot showed protein level of BEX4 in different A498 (F) and PANC-1 (G), respectively. (H) Migration of A498 and PANC-1 cells was assessed by transwell assay after overexpression of BEX4. Scale bar = 100 µm. (I,J) Proliferation of A498 (I) and PANC-1 cells (J) was assessed by cell counting kit-8 (CCK-8) assay after overexpression of BEX4. **p 0.01, ***p 0.001. OE, over expression.