Firstly, the flow chart of this study is shown in Figure 1. Then, the mRNA expression data, copy number alteration threshold data, masked copy number segmentation data, and DNA methylation 450K data of both tumor and normal tissues in the TCGA pan-cancer cohort were obtained from the Firehose database (http://gdac.broadinstitute.org) (Deng et al., 2017). TCPA, mutation, and clinical data were acquired from the UCSC Xena database (https://xenabrowser.net/datapages/) (Wang et al., 2022). The mutation frequency of TUBB in the TCGA cohort was calculated using the cBioPortal database (https://www.cbioportal.org/) (Cerami et al., 2012). Based on the data from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) (Barrett et al., 2013), the TUBB expression at the transcript level was validated. We analyzed and estimated its expression at the protein level based on the CPTAC database (https://proteomics.cancer.gov/programs/cptac) (Zhang et al., 2014). Also, various immune infiltrating algorithms from the TIMER2.0 database (http://timer.cistrome.org) (Li et al., 2020) were adopted to depict the correlations between TUBB expression and tumor immune microenvironment (TIME). Via 99 single-cell datasets from the TISCH database (http://tisch.comp-genomics.org/home) (Sun et al., 2021), immune infiltration results by whole transcriptome analysis were verified. From the CellMiner (https://discover.nci.nih.gov/cellminer/CellMiner) (Shankavaram et al., 2009), GDSC (https://www.cancerrxgene.org/) (Yang et al., 2013), and CTRP (http://portals.broadinstitute.org/ctrp/) databases (Rees et al., 2016), we gained relevant chemotherapy data to illustrate how TUBB expressions interact with drug sensitivity. It should be noted that these public databases are free and open. The study strictly follows the data extraction policy of the databases and does not require the ethical review and approval of the ethics committee.

To observe whether imbalanced expressions of TUBB exist between tumor and normal tissues, differential analyses were conducted in terms of three dimensions.

First of all, we combined the data from the Genotype-Tissue Expression (GTEx) project because the TCGA-data was inadequate, to enhance the confidence by expanding the sample sizes. The “wilcox” test was applied to detect the variability (p < 0.05 was considered significant). The “gganatogram” R package was adopted to visualize the expressions in human’s different organs. Then, the “wilcox” analysis was performed to compare the differential mRNA expressions of TUBB between tumor and normal tissues based on TCGA data. Moreover, we performed the “wilcox” test on paired samples from the TCGA cancer groups (p < 0.05 was considered significant, *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001). Using the “pROC” R package, the importance of TUBB in the diagnosis of pan-cancers was addressed. The area under curve (AUC) values were calculated (between 0.5 and 1, the closer the AUC value came to one, the better diagnostic performance TUBB would have). Usually, AUC values of 0.5–0.7, 0.7–0.9, and >0.9 indicate low accuracy, certainty accuracy, and high accuracy, respectively.

Based on the GEO database, we conducted the external validation at the transcriptional level. At the protein level, the protein expressions of TUBB were validated via the CPTAC database. What’s more, the immunohistochemical (IHC) staining of tissues was collected from the HPA database. The “Kruskal” test was used to identify the expressions of TUBB in different stages or molecular subtypes.

Survival data was obtained from the TCGA database, the “survival”and “survminer” R packages were utilized to describe the relations between the expression of TUBB and prognostic indexes of TUBB (including overall survival (OS), disease-specific survival (DSS), progression-free interval (PFI), and disease free interval (DFI)). Integrating the “Kaplan-Meier (KM)” analysis and the univariate COX analysis, we evaluated whether the TUBB was a protective or risky factor, and finally, the survival landscape of TUBB was drawn with a high level of confidence. Noticeably, the optimal cut-off between the high- and low-TUBB mRNA expression groups was determined via the “survminer” R package when the KM analysis was conducted. The “survfit” function was applied to conduct the “log rank” test to estimate the significance of the high- and low-TUBB expression groups. At last, the “forestplot” R package was used to visualize the results of survival data by the COX analysis.

The cBioPortal (http://www.cbioportal.org) website is a powerful tool for retrieving, downloading, analyzing, and visualizing cancer genomic data from various types of genomic data such as somatic mutations, DNA copy-number alterations (CNAs), and DNA methylation. It can undertake multiple analyses, including mutation analysis and its visualization (Cerami et al., 2012). On cBioPortal, “TCGA Pan Cancer Atlas Study” was selected in the “Quick Search” section (Gao et al., 2013), and “TUBB” was typed to search for its genetic altering traits. The mutation frequency, types, and CNA data were collected from the “Cancer Type Summary” part. The mutation locations of TUBB were shown in a two-dimensional (2D) diagram of the protein structure using the “mutation” function. Somatic copy number alterations (SCNAs) and mutations can increase the CNA of the gene by amplifying and deleting the heterozygosity and homozygosity. Generally, high-frequency SCNA is defined with >5% mutation frequency. Next, the “Spearman” correlations between the expression level of TUBB and CNA scores were calculated to estimate the relation between SCNA and TUBB expression. The “IlluminaHumanMmethylation-450kanno.ilmn12.hg19” R package from “Bioconductor” was conducted to annotate the methylation probe of the TUBB promoter. Through the “Wilcoxon” rank test was performed to detect the differential methylation of TUBB between tumor and normal tissues. A p-value cutoff of 0.05 was used to identify cancers that were significantly hypomethylated or hypermethylated. The Spearman correlation between TUBB expressions and promoter DNA methylation Beta was calculated (p-value < 0.05 was considered significant).

To understand TUBB-related pathways, we divided tumor samples of each type according to TUBB expressions (top 30% and bottom 30%). Then the Gene Set Enrichment Analysis (GSEA) was carried out to compare differential activation or inhibition conditions of 50 hallmark gene sets and 83 metabolic gene sets in different tumors between the high- and low-TUBB group. Yuan H et al. have redefined 14 functional states of malignant tumoral features (Yuan et al., 2019). And, the “z-score” algorithm was proposed by Lee et al. (2008). It can reflect the activity of a given pathway by integrating expressions of characteristic genes. Using the “GSVA” R package, we conducted the “z-score” algorithm on the 14 functional state gene sets. The values were set as the z-score of each gene set. Then, the “Pearson” correlation analysis was performed to calculate statistical relations between TUBB expression and the z-score of each gene set. Additionally, we identified the differential genes between the high- and low-TUBB group. Moreover, to search for transcription factors that may affect the expression of TUBB, we used the CistromeDB database (http://dbtoolkit.cistrome.org/) to identify potential regulatory upstream factors for TUBB. As known, protein-protein interacting data usually contains those unlikely biological interactions that are impossible to happen in living cells. Therefore, we used the ComPPI database (https://comppi.linkgroup.hu/) which introduces two novel quantitative scores, the Localization Score and the Interaction Score, describing the calculated probability of the data correctness, to gain the proteins that interact with TUBB. TUBB has been identified as a risky factor according to four survival indexes in sarcoma (SARC). Hence, we conducted the GSEA on multiple gene sets to analyze possible pathways that SARC may participate in. It was believed that if one gene had differential expression in more than five cancers, its functional gene would be thought to be related to TUBB and undertake a KEGG enrichment analysis to recognize the conservative functions or pathways that TUBB takes part in pan-cancers. Subsequently, we adopted the permutation test to identify the mutation that was statistically related to TUBB expressions in SARC. Then, using Fischer’s exact test, we estimated the relationship between E2F1 and TUBB expressions. In the end, we conducted the Spearman correlation analysis to systemically identify BEND3-related proteins.

The continuous interactions between tumor cells and the tumor microenvironment (TME) play a decisive role in tumor initiation, progression, metastasis, and responses to therapies. As a consequence, we confirmed correlations between TUBB expressions and immune-related genes (immunostimulatory genes, immunoinhibitory genes, chemokine genes, chemokine receptor genes, and MHC genes). We also adopted the same approach to conduct the differential analysis on TIP scores. Next, from the TIMER2.0 databases, seven most advanced algorithms were applied (CIBERSORT, CIBERSORT_ABS EPIC, MCPCOUNTER, QUANTISEQ, TIMER, XCELL) to estimate the immune filtration profiles of TCGA cancers (Li et al., 2020). At last, we used the TISCH database to download the expression landscape of TUBB in 99 single-cell datasets of 38 tumors, which helped us verify the TME analytic result at the single-cell level. In summary, we provided a comprehensive analysis and visualized the TME and immune infiltration profiles of TUBB in a pan-cancer cohort. It is worth noting that, we analyzed the anti-tumor immunity in seven steps of the cancer-immunity cycle (release of cancer cell antigens, cancer antigen presentation, priming and activation, trafficking of T cells to tumors, infiltration of T cells into tumors, recognition of cancer cells by T cells, killing of cancer cells (Chen and Mellman, 2013) and conducted the Gene Set Variation Analysis (GSVA) on each step, and then we compared their differences in the high- and low-TUBB groups.

Based on the GSCA database (http://bioinfo.life.hust.edu.cn/web/GSCALite/), we analyzed the relations between TUBB expressions and drug sensitivity (Liu et al., 2023). GSCA website contains 750 small molecule drugs from GDSC and CTRP databases. Besides, the gene expression data is used to mine the valuable small molecule drugs related to it. In addition, the National Cancer Institute (NCI) database was utilized. NCI serves as a trusted source of cancer information, and the platform of cancer cell lines from NCI has been widely used to screen out drugs that are related to a certain gene’s expression. NCI-60 is a collection of 60 cancer cell lines from nine different kinds of cancers (leukemia, colon, lung, central nervous system, kidney, melanoma, ovarian, breast, and prostate cancers). Its data comes from the “CellMiner” database (Reinhold et al., 2023), for us to analyze the Spearman correlations between TUBB mRNA expressions and the z-score of drug sensitivity. Moreover, differentially expressed genes between the high- and low-TUBB group in different cancers were identified. We collected the top 150 upregulated and downregulated genes as TUBB-related biomarkers. The CMAP_gene_signatures.RData document containing 1,288 compound-related characteristics was downloaded from the BHKLAB database (https://www.pmgenomics.ca/bhklab/sites/default/files/downloads) and used to calculate matching scores. All analyzing processes kept to the methods from the previous literature (Malta et al., 2018). The R software was used to summarize and display the top five results of 32 kinds of cancers.

All data was processed using the web tools and R software (V.4.3.0, Institute of Statistics and Mathematics, Vienna, Austria). The Pearson correlation analysis was used on normally distributed data, otherwise the Spearman correlation analysis was used. The Kruskal–Wallis rank sum test, Wilcoxon rank sum and Signed rank tests were used to detect the differences between multiple variables or two variables, respectively. Using the “survival” R package, the COX and KM survival analyses were conducted. KM method adopted the log-rank test to detect significance. The “survminer” R package was used to visualize results from KM analysis. Hazard ratio (HR) and 95% confidence interval (CI) were used to describe relative risks. The “pROC” R package was utilized to perform ROC analysis to estimate the diagnostic ability of TUBB. All statistical tests were two-tailed. p-value < 0.05 was deemed as statistically significant. p-value < 0.0001 was deemed as greatly statistically significant (*p < 0.05, **p < 0.01, ***p < 0.001, and ****p < 0.0001).

To identify patterns of TUBB regulations in cancers, we combined TCGA and GTEx data to expand sample sizes and gained a boxplot describing TUBB expressions in pan-cancer (Figure 2A). Figure 2B showed TUBB expressions in different organs. TUBB had a differential expression in most cancers. It was greatly upregulated across cancers. Then, the differential analysis based on TCGA samples and paired samples were shown in Figures 2C, D. With the aid of logistics regression analysis using TCGA and TCGA-GTEx data, above results were validated (Figure 2E). An external validation of TUBB mRNA expressions was based on GEO database (Supplementary Figure S1). From CPTAC database, the validation at the protein level was conducted (Figures 2F–K). It was observed that results from different angles and databases had good consistency. IHC results showed that, most cancer cells displayed weak to moderate cytoplasmic immunoreactivity. A few seminomas and carcinoids were strongly stained. Most hepatocellular carcinomas were negative (Figure 2L). Moreover, TUBB expressions were correlated with tumor stages in seven cancers (ACC, ESCA, KIRC, KIRP, LIHC, SKCM, STAD) (Supplementary Figure S2), implying that TUBB may have a relationship with the progression of some cancers. The estimated ROC curves showed that (Supplementary Figure S3), the TUBB mRNA expression level showed satisfactory sensitivity and specificity for the diagnosis of nine kinds of tumors (AUC>0.7). After expanding normal sample sizes with GTEx data, the results also were robust (Supplementary Figure S4). The repeatable and consistent results were demonstrated across multiple databases, multiple tumors, multiple methods, and multiple omics, suggesting that dysregulation of TUBB expression may play a crucial role in different cancers and was highly unlikely to be a false finding due to technical artifacts, chance, or bias in the eligibility criteria for TCGA samples. In addition, by the chi-square test, we found that C1 and C2 immune subtypes were more in the high-TUBB group (Figure 3A). The same results were obtained using the Kruskal test, that is, TUBB was expressed more in the C1 and C2 immune subtypes compared to other subtypes (Figure 3B). Interestingly, TUBB showed differences in a large number of molecular subtypes. For example, it had the lowest expression in LumA group, and in the basal group, its expression was the highest in breast cancer (BRCA), meaning that it had great values in the precision molecular stratification therapies and prognosis prediction (Figure 3C).

To further dissect the clinical relevance of TUBB in cancers, the role of TUBB was analyzed. Survival profiles of pan-cancers showed that TUBB expressions were related to many kinds of survival indicators in pan-cancers (Figure 4A), and the relations had homogeneity. TUBB was identified as a risky factor in cancers, especially in KIRP, PAAD, and SARC because TUBB was risky weighed by all indicators DFI, DSS, OS, and PFI. In several cancers, TUBB was protective. For example, UVM patients with higher TUBB expressions would have better survival probability. Due to the highly heterogeneity across various cancers, in some cancers (BLCA, CHOL, ESCA, UCS), the tendency was vacant. Therefore, TUBB may play different roles in pan-cancers, suggesting further explorations of it should be addressed. Forest plots showed COX survival analysis results of four survival indicators. The hazard ratios of each cancer were listed, too (Figures 4B–E). KM curves were used to display the log-rank test results in KIRP, PAAD, and SARC (Figure 4F).

We analyzed genomic data (genetic variation, somatic copy number alteration (SCNA), mRNA expression, and DNA methylation) of tumor and normal tissues from the TCGA cohort. Mutation sites of TUBB were visualized in 2D and 3D graphics (Figures 5A, B). TUBB had genetic changes in most cancers, and the most common alteration types were amplification and mutation (Figure 5C). Obviously, SCNA plays an essential role in regulating gene expression in cancers (Figure 5D). To study genetic changes of TUBB in cancers, we checked the percentage of SCNA. SCNA appeared frequently in most cancers (more than 5% among all samples). Only in a few types, the frequency was low (Figure 5E). Therefore, we estimated how SCNA affected TUBB mRNA expressions by calculating the Spearman correlation between TUBB expressions and masked copy number segment of TCGA. In most cancers, the mRNA expression of TUBB was significantly correlated with SNCA (Figure 5F). Therefore, it could be concluded that TUBB CNAs were common among cancers, and it could affect TUBB expressions. Besides, promoter methylation of TUBB can also regulate gene expressions. And, abnormal DNA methylation of promoters was associated with tumorigenesis. The methylation patterns of TUBB in pan-cancers were consistent with each other. From all results, it was observed that TUBB mRNA expressions were negatively correlated with DNA methylation (Figure 5G), and higher methylations existed in tumor tissues more than in normal tissues in most cancers (Figure 5H) because when TUBB is significantly overexpressed, the body increases promoter methylation to overcome this dysregulation to keep the balance. Great differences in the methylation patterns of TUBB suggested the complexity of TUBB regulation and the specificity of this process among different cancer species. In addition, we selected 10 transcription factors with the highest scores to annotate peaks and found that they were usually located in the promoter region of TUBB, showing that they may regulate the expression of TUBB.

Samples with top 30% and bottom 30% TUBB expressions were set as the high- and low-TUBB group. Based on two groups, we conducted the GSEA to explore relevant cell signaling in cancers. Furthermore, metabolism-related pathways were analyzed systematically, and the results had satisfying consistency across cancers, indicating the functions of TUBB were highly conservative. What’s more, pathways related to cell-cycle were enriched in tumors with higher TUBB. Besides, TUBB may take part in metabolic disorder processes in cancers. Meanwhile, TUBB was believed to inhibit the course of drug metabolism, thus it may be related to responses to chemotherapy (Figure 6A). We analyzed the correlations between z-scores of 14 symbolic functional states of cancers (angiogenesis, apoptosis, cell cycle, differentiation, DNA damage, DNA repair, EMT, hypoxia, inflammation, invasion, metastasis, proliferation, quiescence, stemness) with the z-score of TUBB expressions. Among them, the z-score of the cell cycle had the highest R-value (0.51), showing a positive relation with TUBB expressions (Figure 6B). Then, in multiple cancers, we ran a Pearson correlation analysis between the z-score of TUBB expressions and the z-score of the cell cycle state. Learning from Figure 5C, in thymoma (THYM) the highest coefficient of 0.87 was detected (Figure 6C). CHIP-seq demonstrated that E2F1 might be an upstream transcription factor modulating the TUBB expression (Figure 6D). ComPPI helped identify the genes that may interact with TUBB (Figure 7A). The KEGG enrichment analysis was conducted on genes that were highly expressed in the high-TUBB group. Figure 7B showed that TUBB had relations with a lot of functions, especially with signal transduction, immune, and cancer-related pathways. As previously mentioned, TUBB was determined as a risky factor in terms of four survival indexes in SARC. Then we performed the GSEA on various gene sets to fully analyze the pathways that TUBB may be involved in (Figure 7C). From Figure 7D, the expression of TUBB probably was relevant to ATRX mutations. In SARC, E2F1 was statistically associated with TUBB expressions (Figure 7E). Besides, based on the TCPA database, some TUBB-related proteins were identified (Figures 7F–O). Among them, CYCLINB1 had the most significant relation with TUBB mRNA expressions with the highest coefficient of 0.57.

First of all, the heatmap displayed the landscape of immune infiltration of TUBB across cancers (Figure 8A). TUBB was significantly negatively correlated with immune-related genes (MHC, immune-inhibitor, immune-stimulator, chemokines). Notably, we scored seven steps and compared the differences between the high- and low-TUBB group, and the tendency was consistent in pan-cancers, namely the TIP scores were lower (Supplementary Figure S5). To illustrate cells regulated by TUBB in the TME, we used the TIMER2.0 database to explore correlations between TUBB mRNA expressions and immune infiltration and stroma cell abundance (Figure 8B). For example, in almost every cancer the TUBB mRNA expression level was positively related to the abundance of CD4⁺ Th1 cells. In SARC, the TUBB mRNA expression level was negatively related to the abundance of B cells and CD8⁺ T cells across different software, which could possibly explain why TUBB was risky in cancers. Different proportions of immune infiltration and unique TMEs in pan-cancers existed, so the correlations showed various change laws. However, we utilized seven different algorithms to conduct the analysis and made the results verify each other to ensure the accuracy of the study. Additionally, we could see from Figure 8C that the single-cell analysis showed that though TUBB was not expressed strikingly in most tumors, it mainly originated from malignant cells and proliferative T cells. This echoed previous results. In other words, TUBB might encourage the formation of immunological rejection or “immunological desert”, and play a vital role in immunity-cancer crosstalk, especially in immune escape.

Using Cellminer data, we found that TUBB was positively related to the sensitivity of a lot of drugs (Figure 9A), but negatively related to the half maximal inhibitory concentration (IC50) values of many drugs based on CTRP, GDSC databases (Figures 9B, C). Therefore, we could deduce that TUBB was a potential chemotherapy-sensitive gene. To further explore underlying therapeutic regimens that offset the tumorigenic effect mediating by TUBB, CMap analysis was employed. Thus, a TUBB-related signature containing 150 significantly upregulated and 150 significantly downregulated genes that were selected by screening patients with high- or low-TUBB expressions in various cancers, was generated. Using the Extreme Summarization (XSum) method, TUBB-related traits were compared with CMap gene traits, and the similarity scores for 1,288 compounds were obtained. Fasudil, W.13, AH.6809 and X4.5. dianilinophthalimide showed relatively lower scores in most cancers, which demonstrated that they perhaps inhibit TUBB-mediated carcinogenic effects (Figure 9D). These findings provided substantial support for the validity of our predictions, although further researches are needed to elucidate the underlying mechanisms.

To verify the expression of TUBB among the tissue samples of cancer patients, IHC staining was conducted. Figure 10A–C displayed the negative results of TUBB expression in the control subjects. TUBB exhibited strong diffuse staining in osteosarcoma tissues (Figures 10D, E). Positive staining was also found in chondrogenic sarcoma tissues (Figure 10F).