TMUB1 expression is associated with the prognosis of colon cancer and immune cell infiltration

Acquisition of bioinformatics analysis data

The TCGA database project (https://portal.gdc.cancer.gov/) was used to download colon cancer (colon) identified by high-throughput sequencing RNA data (each 1 million base fragments (fragments per kilobase million (FPKM) format) and the corresponding clinical pathological information data. Herein, RNA sequencing data was converted from the FPKM format to the transcripts per million reads format. The study was conducted in accordance with TCGA recommendations, and informed consent was obtained before data collection. The Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/) of the National Centre for Biotechnology Information was used in this study. A colon cancer sample group and a control group (healthy tissue, adjacent non-cancer tissue), with at least 10 samples in each group, were included in the study. Finally, four gene expression profiles (GSE83889, GSE9348, GSE23878, and GSE47756) were selected. R (version 3.6.3) was used for data normalisation and statistical analysis and visualisation. The HPA database (https://www.proteinatlas.org/) retrieved TMUB1 genes for normal colorectal tissue and colorectal cancer tissue. The Institutional Review Board of Shunde Hospital Southern Medical University (The First People’s Hospital of Shunde), approved this study with the approval number 20210730.

Clinical significance of TMUB1 expression in colon cancer

Receiver operating characteristic (ROC) analysis was used to compare TMUB1 expression in colon cancer and the adjacent tissues, and the predictive value of TMUB1 in colon cancer diagnosis was tested. The clinical prognostic information of patients with colon cancer included overall survival (OS) and disease-specific survival (DSS). The prognosis was analysed using Kaplan–Meier (K–M) analysis and univariate and multivariate Cox regression analysis. R package: Survival (version 3.2.10) was used for statistical analysis of survival data (Liu et al., 2018). R package (rms) was used to construct, analyse, and visualise lipograms and calibration diagrams.

Integration of differentially expressed gene screening and construction of a protein–protein interaction (PPI) network

RNA sequencing data in level 3 HTSeq-Counts format for differentially expressed gene (DEG) was identified using R package DESeq2 (version 1.26.0) and ggploT2 (version 3.3.3) (Love, Huber & Anders, 2014). A protein–protein interaction (PPI) network was constructed by analysing the molecular interactions in the TMUB1 protein list through the STRING database (https://cn.string-db.org/), and the relationship of the protein network was displayed (Szklarczyk et al., 2019). R package (igraph) was used for visualisation (Mora & Donaldson, 2011).

TMUB1 functional annotation and Kyoto Encyclopaedia of Genes and Genomes pathway analysis

Gene Ontology (GO) functional annotation and Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed for TMUB1. Statistical analysis and visualisation were performed using version R3.6.3. Enrichment analysis was performed using clusterProfiler (version 3.14.3), and visualisation was performed using ggploT2 (Yu et al., 2012).

TMUB1 differential GSEA

Herein, R package’s clusterProfiler was used to perform GSEA to clarify significant functional and pathway differences between the high TMUB1 and low TMUB1 groups. Data sets from MSigDB after the adjustment of P < 0.05 error detection rate, false discovery rate of P < 0.25, and enrichment of standardised scores (—NES—) > 1 were considered significantly enriched.

Correlation analysis of TMUB1 immune infiltration

Immune infiltration analysis of TMUB1 was performed using the single-sample GSEA algorithm in R package, Gene Set Variation Analysis (Bindea et al., 2013). Twenty-four types of infiltrating immune cells were obtained. The correlation between TMUB1 and the enrichment scores of the 24 types of immune cells was analysed using Spearman’s correction method. The enrichment scores of TMUB1 high and low expression groups were analysed using the Wilcoxon rank-sum test.

Cell culture

Human intestinal epithelial cells (NCM460) were cultured in Roswell Park Memorial Institute 1640, comprising 10% foetal bovine serum. Human intestinal carcinoma cells (SW480, RKO, HCT116, LoVo, and HT29) were cultured in Dulbecco’s modified Eagle medium with 10% foetal bovine serum and were incubated at 37 °C, under 5% CO₂ and 95% humidity. All cells were procured from CAS Cell Bank.

Tissue specimens

Cancer tissue and adjacent tissue samples were collected from 10 cases. All clinical samples were stored in a liquid nitrogen tank after quick-freezing within 40 min in vitro. All tissue samples used for total RNA extraction had RNA protectants added to them. All tumor samples were approved by the ethics committee of our university, and informed consent was obtained from all patients and received written informed consent from participants of the study. A clinicopathological diagnosis of colorectal cancer was confirmed in the included specimens.

RNA extraction

Total RNA was extracted from tissues using TRIzol reagent (TaKaRa, Tokyo, Japan) as per the manufacturer’s instructions. PrimeScript RT reverse transcription kit (#RR047A, TaKaRa, Tokyo, Japan) was used to reverse transcribe the extracted RNA. The amount of isolated RNA was detected using the SYBR green polymerase chain reaction (PCR) kit (#RR820A, TaKaRa, Tokyo, Japan). Quantitative reverse transcription (qRT)-PCR was performed using the cfX-96 coupled real-time fluorescence quantitative system (No. 788BR07388 Bio-Rad, Hercules, CA, USA).

The primer sequences used were as follows:

5′-CTACCTCATGAAGATCCTCACCGA-3′(beta-actin, forward)

5′-TTCTCCTTAATGTCACGCACGATT’ (beta-actin, reverse)

5′-GTGTCCACGAGAGTCGGTC’ (TMUB1, forward)

5′-AGGGCCGGTACTGGATCTG-3′(TMUB1, reverse)

Statistical methods

Statistical analysis was performed using GraphPad Prism 8. The two sets of data were compared using an unpaired two-tail test. The chi-square test and rank-sum test were used to compare parameters. The correlation between variables was evaluated using Spearman’s correlation analysis. The relationship between the TMUB1 expression level and clinical features was analysed using univariate logistic regression. P < 0.05 was considered statistically significant (*, P < 0.05; **, P < 0.01; ***, P < 0.001).

TMUB1 expression in colon cancer and other tumors

The TCGA database showed that TMUB1 mRNA levels were highly expressed in different tumors. Among the 33 tumor types, TMUB1 was significantly overexpressed in 23 tumors (Fig. 1A), particularly in gastrointestinal tumors (cholangiocarcinoma, colon adenocarcinoma, oesophageal carcinoma, liver hepatocellular carcinoma, rectum adenocarcinoma, and stomach adenocarcinoma). Similarly, TMUB1 was found to be significantly overexpressed in gastrointestinal tumors in paired samples (Fig. 1B). Particularly, TMUB1 expression was significantly higher in colon cancer tumors than in pericancerous tissues (Fig. 1C) and in paired tumors than in normal tissues (N = 41) (Fig. 1D).

TMUB1 expression in the HPA and GEO databases

The HPA database revealed that TMUB1 was significantly overexpressed in colon cancer tissues. In the HPA database, deeper immunohistochemical staining of cancer tissues compared with normal tissues was observed, suggesting that TMUB1 was significantly overexpressed in colon cancer tissues (Fig. 2A). The GEO database data sets (GSE9348, GSE23878, GSE47756, and GSE83889) revealed that TMUB1 expression in colon cancer was significantly higher than that in normal tissues (Figs. 2B–2E).

Clinical characteristics of TMUB1 in COAD and prediction of its diagnostic and prognostic value

We studied the clinicopathological features of TMUB1 expression in COAD patients with differential topoisomerase II β-binding protein 1 interacting checkpoint and replication regulator expression, as shown in Table 1. The value of TMUB1 in the differential diagnosis of COAD was proven using the ROC curve. An area under the curve (AUC) of 0.845 indicated that TMUB1 was highly sensitive and specific in diagnosing colon cancer (Fig. 3A). K–M analysis verified the prediction of clinical survival using TMUB1, with a statistical significance in OS (Fig. 3B) (hazard ratio (HR) =1.98, P = 0.002) and DSS (Fig. 3C) (HR =2.23, P = 0.003). A poor prognosis was observed in colon cancer patients with a high TMUB1 expression than in those with a low TMUB1 expression. Table 2 presents the correlation between TMUB1 expression and the clinicopathological features analysed using logistic regression. High TMUB1 expression was associated with a high risk ratio in M-stage COAD (odds ratio (OR) = 2.082 [1.210–3.675], P = 0.009), clinicopathological advanced stage (OR = 1.487 [1.029–2.153], P = 0.035), and lymphatic metastasis.

The significance of TMUB1 in the clinical prognosis using univariate Cox regression analysis

The results of the univariate Cox regression analysis results based on the OS (Fig. 4A) and Table S1 showed that high TMUB1 expression was associated with a higher risk of clinically adverse prognostic factors. TMUB1 was a significant risk factor for clinical T3 and T4 stages (HR = 3.072, P = 0.004); N1 and N2 (HR = 2.592, P < 0.001); M1 (HR = 4.193, P < 0.001); Stages III and IV (HR = 2.947, P < 0.001); the presence of lymphatic invasion (HR = 2.450, P < 0.001). The prognosis based on the DSS has the same results (Fig. 4B). Further multivariate Cox regression analysis of the OS showed (Fig. 4C) revealed that high TMUB1 expression was an independent risk factor for OS (HR = 1.917, P = 0.009). Similarly, high TMUB1 expression was observed in clinical T3 and T4 stages (HR = 3.332, P = 0.047), M1 (HR = 1.906, P = 0.03), and stages II and IV (HR =4.232, P = 0.021). Another independent risk factor was the presence of lymphatic invasion (HR = 1.771, P = 0.034). TMUB1 had a significant predictive power for DSS (Fig. 4D, Table S2).

The predictive value of TMUB1 in the diagnosis and prognosis of colon cancer

Clinical T stage, M stage, pathological stage, lymphatic metastasis, and TMUB1 expression were constructed, and calibration curves were drawn to test the validity of the OS rolograms. The results revealed that the c-index of the OS rate of colon cancer was 0.747 (Fig. 5A). This helped predict the OS rate of colon cancer at 1 and 3 years in addition to predicting the prediction OS rate at 5 years (Figs. 5B, 5C, and 5D). Similarly, the M stage, the pathological stage, lymphatic metastasis, and TMUB1 expression were used to construct the linear DSS prediction graph, and the c-index was 0.813 (Fig. 5E). Calibration curves revealed relatively accurate predictions of the clinical outcomes at 1, 3, and 5 years (Figs. 5F, 5G, and 5H). Similarly, in Fig. S1, DSS is also well predicted (Fig. S1).

TMUB1 in the clinical prognosis subgroup of colon cancer

The prognostic subgroups of TMUB1 in terms of colon cancer OS and DSS were further evaluated. In the clinical subgroup K–M OS analysis, clinical T3 and T4 staging (HR = 2.35, P = 0.001), pathological stage (HR = 1.95, P = 0.010), age (HR = 2.08, P = 0.003), body mass index (BMI) (HR = 3.21, P = 0.011), and body weight (HR = 8.71, P = 0.041) revealed a poor prognosis of TMUB1 with high expression (Fig. 6A). Similarly, in the clinical subgroup K–M DSS analysis, clinical T3 and T4 staging (HR = 2.61, P = 0.001), N1 and N2 (HR = 2.06, P = 0.018), pathological staging (HR = 2.49, P = 0.003), and age (HR = 2.27, P = 0.015) revealed a prognosis of TMUB1 with high expression (Fig. 6B).

DEGs, co-expressed genes, and a PPI network of TMUB1 in the TCGA in colon cancer

Variance analysis with volcanic figure illustrated —log2 (FC)— > 1 and p.adj < 0.05. Number of molecules at the threshold of 0.05. Volcano shows meet —log2 (FC)— > 1 & p.adj < 0.05. The 0.05 threshold had 1172 ids. Under this threshold, 128 ids had high expression (logFC positive) and 1,044 ids had low expression (logFC negative) (Fig. 7A). The top 10 genes positively and negatively associated with TMUB1 are illustrated in the co-expression heat map (Fig. 7B). The top 20 interacting proteins associated with TMUB1 by constructing PPI protein interaction networks are as follows: transmembrane protein (TMEM) 126A, autocrine motility factor receptor, family with sequence similarity 3 metabolism-regulating signalling molecule C, TMEM209, endoplasmic reticulum lipid raft-associated protein (ERLIN) 2, glutamate receptor interacting protein 1, saccharopine dehydrogenase, calcium-modulating ligand, ERLIN1, ring finger protein 139, TMEM53, protein interacting with C kinase-1, cleft lip and palate transmembrane protein 1, ubiquitin-associated domain-containing protein 2, fas associated factor family member 2, calsyntenin 3, E3 ubiquitin-protein ligase synoviolin, chromosome 17 opening reading frame 62, NCLN, and nurium (Fig. 7C).

Functional enrichment analysis of TMUB1

By identifying differential genes, functional annotation analyses were performed, including bioprocess, cellular component, and molecular function of three types of GO and the KEGG pathway enrichment analyses. The findings of the GO analysis revealed that TMUB1 was involved in ion transport, ion channel regulation, intercellular adhesion, and extracellular matrix composition (Fig. 8A). The KEGG analysis revealed that the neuroactive ligand–receptor interaction pathway, the calcium signalling pathway, the cyclic adenosine monophosphate signalling pathway, and the extracellular matrix (ECM)-receptor signalling pathways such as the interaction pathway were significantly enriched (Fig. 8B). The following TMUB1-related signalling pathways were identified by GSEA through comparative differential gene analysis: the insulin pathway, bladder cancer, glutathione metabolism, deoxyribonucleic acid (DNA) methylation, signalling by notch, and cell cycle checkpoints (Fig. 8C).

Correlation analysis of TMUB1 immune function in colon cancer

The degree of TMUB1 immune infiltration in colon cancer was examined, and it was discovered that TMUB1 was associated with a variety of immune cells. Figure 9A shows that TMUB1 was positively correlated with CD56^bright natural killer cells and other cells (R = 0.268, P < 0.001), and negatively correlated with T helper cells (R =  − 0.407, P < 0.001), central memory T (Tcm) cells (R =  − 0.387, P < 0.001), and macrophages (R = −0.238, P < 0.001). The immune infiltration and matrix scores were evaluated after TMUB1 was divided into high-low expression groups, and the results revealed that the immune infiltration and matrix component scores of high-expression TMUB1 were lower (Fig. 9B). Meanwhile, high-expression TMUB1 in colon cancer showed lower infiltration levels in T helper cells, Tcm cells, T helper 2 (Th2) cells, and macrophages (Fig. 9C).

TMUB1 methylation level in colon cancer

We assessed the TMUB1 methylation levels in colon cancer. The results revealed lower levels of methylation in colon cancer compared with normal tissues (Fig. 10A). Through Illumina methylation 450 (high-throughput methylation chip detection platform) and Spearman’s correlation analysis, it was observed that in RNA sequencing data of colon cancer, TMUB1 expression was negatively correlated with the beta corresponding to the methylation probe CG7847101 (R =  − 0.329, P < 0.001) (Fig. 10B). Among other features, lower TMUB1 methylation levels in colon cancer were observed in TP53 mutations, adenocarcinoma categories, clinicopathological stage III/VI, and increased age (Figs. 10C–10F).

qRT-PCR for verifying TMUB1 expression levels in different colorectal cancer cell lines and tissues

To determine their expression levels in comparison to normal intestinal epithelial cells, human intestinal epithelial cells (NCM460) and human intestinal carcinoma cells (SW480, HCT116, RKO, LoVo, and HT29) were cultured. It was found that TMUB1 expression in most intestinal carcinoma cells was higher than that in normal NCM460 cells (Fig. 11A). At the same time, paired colorectal cancer tissues from 10 patients were collected. Of the 10 paired colorectal cancer tissues (Fig. 11B), eight pairs had TMUB1 expression levels higher than normal, while the remaining two pairs showed no difference (NATs: Normal Tissue Around tissues; tumor tissue).