Mito-fission gene prognostic model for colorectal cancer

Data extraction

RNA-seq, clinical, and survival data were retrieved from The Cancer Genome Atlas (TCGA) database, encompassing 606 CRC tissue samples (including survival data) and 51 normal tissue samples, hereafter referred to as TCGA-CRC. Additionally, the GSE103479 dataset was extracted from the Gene Expression Omnibus (GEO) database, which included 155 CRC tissue samples with associated survival data. Two gene sets, GOBP MITOCHONDRIAL FISSION and GOBP POSITIVE REGULATION OF MITOCHONDRIAL FISSION, were selected from MSigDB using “mitochondrial fission” as the keyword, resulting in the identification of 40 mitochondrial fission-related genes (MFRGs) for subsequent analysis (Table 1).

Identification of DEGs

In TCGA-CRC samples, differentially expressed genes (DEGs) were identified using the DESeq2 package (version 1.34.0) (Love, Huber & Anders, 2014), with thresholds set at P.adj < 0.05 and log₂FC > 1.5. The top 20 DEGs, sorted by log₂FC, were visualized via volcano and heat maps.

Analysis of MFRGs scores

MFRG scores for each sample were calculated using the GSVA package (version 1.42.0) (Hänzelmann, Castelo & Guinney, 2013) to assess differences between sample groups. A rank-sum test was employed to evaluate differences in MFRG scores between the two sample groups. The impact of MFRGs on the survival of patients with CRC was evaluated by stratifying samples based on the optimal cut-off values of MFRG scores. Kaplan–Meier (KM) curves were plotted for high- and low-MFRG score groups using the survminer package (version 0.4.9) (Li et al., 2020). in order to understand the relationship between MFRGs scores and clinical characteristics (age (50-year cut-off), gender, ethnicity and T/N/M stage), CRC patients were divided into different clinical subgroups, and the differences in MFRGs scores between different subgroups were compared using the Wilcoxon test.

Weighted gene co-expression network analysis

The expression matrix of all CRC samples from TCGA-CRC was analyzed using weighted gene co-expression network analysis (WGCNA) (version 1.70.3) (Langfelder & Horvath, 2008). Initially, outlier samples were excluded by clustering. The optimal soft threshold (R² = 0.85) was determined by evaluating the relationship between the soft threshold and the scale-free network evaluation coefficient. A phylogenetic tree of genes was then constructed based on gene similarity. Using MFRG scores as the phenotype, the correlation between the phenotype and gene modules was computed. Correlation coefficients and P-values were calculated, and a correlation heat map was generated. The module with the strongest correlation to the MFRG score (R > 0.5 and P < 0.05) was selected as the key module. Further gene filtering within the module was performed by setting the gene significance (GS) to 0.4 and module membership (MM) to 0.4, thereby identifying key module genes associated with the MFRG score.

Acquisition and enrichment analysis of DE-MFRGs

DE-MFRGs were identified by intersecting DEGs and the key module genes associated with MFRGs using the VennDiagram package (version 1.7.1) (Ito & Murphy, 2013). To explore the biological functions of the DE-MFRGs, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed using the clusterProfiler package (version 4.2.2) (Wu et al., 2021), with a significance threshold set at P < 0.05.

Risk model construction

To assess the prognostic impact of DE-MFRGs on patients with CRC in TCGA-CRC, a risk model was constructed. First, univariate Cox regression analysis was performed using the survival package (version 3.5-3) (Deng & Thompson, 2023) to calculate risk scores and identify significant genes (P-value < 0.05). In the second step, Least Absolute Shrinkage and Selection Operator (LASSO) analysis was conducted using the glmnet package (version 4.1-4) (Friedman, Hastie & Tibshirani, 2010) to select genes with strong correlations for constructing the optimal risk model. In the third step, the proportional hazards (PH) assumption was tested for the genes selected from the LASSO analysisd, and those with P > 0.05 were included in further analysis, followed by multivariate Cox regression analysis. Finally, a stepwise regression method was applied to optimize the multivariate Cox regression results by entering variables one by one and checking their significance, removing non-significant variables, thereby constructing the optimal regression equation and identifying the hub genes associated with prognosis.

Validation of risk model

Using the expression levels of hub genes and the risk coefficients obtained from stepwise regression analysis, the risk score formula ($\text{Riskscore}={\sum }_{\mathrm{i}=1}^{\mathrm{n}}\text{coef}\left(\mathrm{gen}{\mathrm{e}}_{\mathrm{i}}\right)\times \text{expr}\left(\mathrm{gen}{\mathrm{e}}_{\mathrm{i}}\right)$) was used to calculate the risk score for each patient, with patients with CRC categorized into high- and low-risk groups based on the median risk score. The prognostic value of the risk model was assessed by comparing survival statuses between the two risk groups, and KM survival curves were plotted using the survminer package (version 0.4.9) (Li et al., 2020). The model’s validity was further assessed by plotting receiver operating characteristic (ROC) curves. The risk model was then validated in the GSE103479 dataset using the same methodology.

Independent prognostic analysis and clinical characteristic correlation analysis

For univariate Cox analysis, six clinical characteristics (age, gender, race, T/N/M stages, and risk score) were included in the model. The results underwent PH assumption testing and multivariate Cox analysis to identify independent prognostic factors associated with CRC. Additionally, to evaluate the risk model’s applicability for patients with CRC, a nomogram predicting 1-, 2-, and 3-year survival was constructed using the rms package (version 6.5-1) (Sachs, 2017). The nomogram’s validity was verified by calibration curves and ROC analysis. Finally, in TCGA-CRC, differences in risk scores across the six clinical characteristics were analyzed using the rank-sum test (P < 0.05).

GSEA analysis

Differential analysis between the two risk groups in TCGA-CRC was performed using DESeq2 (version 1.34.0) (Love, Huber & Anders, 2014), with —log₂FC— calculated and ranked. GSEA was conducted using the KEGG background dataset (c2.cp.kegg.v2023.1.Hs.symbols.gmt), with a threshold set at FDR < 0.05. The top six enriched signaling pathways were visualized using the enrichplot package (version 1.18.3) (Zhang et al., 2019).

Immune related analysis

To further investigate the differences in the immune microenvironment between the high- and low-risk groups in CRC samples, the relative abundance of 22 immune cell types across all TCGA-CRC samples was evaluated using the “CIBERSORT” algorithm, and immune cell proportions were compared between high- and low-risk group samples. Differences in immune cell composition between the risk groups were assessed with the “Wilcoxon” test (P < 0.05). The expression levels of 48 immune checkpoints were also compared between these groups, and the TIDE score, immune exclusion score, and immune dysfunction score were calculated using the TIDE database (http://tide.dfci.harvard.edu/) to evaluate immune therapy efficacy.

Drug sensitivity analysis

The half maximal inhibitory concentration (IC50) for 138 common chemotherapeutic agents was calculated for all TCGA-CRC individuals using the “pRRophetic” algorithm (version 0.5) (Geeleher, Cox & Huang, 2014). Differences in IC50 values between the two risk groups were compared using the Wilcoxon test. A Spearman correlation analysis was performed to assess the relationship between the sensitivity to chemotherapeutic agents and hub gene expression (|r| > 0.4).

Pan-cancer analysis

To examine the expression of hub genes across different cancer types, the Wilcoxon test was conducted using the rstatix package (version 0.7.2) with a significance threshold of P < 0.05. Results were visualized with the ggplot2 package (version 3.5.1). Box plots displayed differential expression of hub genes in pan-cancer samples, while violin plots illustrated the differences in hub gene expression between tumor and normal samples in TCGA-CRC.

RNA isolation and quantitative real-time polymerase chain reaction (qRT-PCR)

Twenty tissue samples, including 10 CRC tumor and 10 normal samples, were collected, frozen, and stored at −80 °C. Each 50 mg sample was lysed with TRIzol reagent for total RNA isolation, following the manufacturer’s instructions. RNA concentration was measured by NanoPhotometer N50, and its quality was assessed based on concentration, purity (A260/A280 ratio), and the amplification curve. RNA was reverse transcribed into cDNA using the SureScript First Strand cDNA Synthesis Kit (Servicebio, Wuhan, China). The qRT-PCR reactions consisted of three µL of reverse transcription product, five µL of 2xUniversal Blue SYBR Green qPCR Master Mix, and one µL each of forward and reverse primers. The amplification protocol included an initial denaturation step at 95 °C for 1 min, followed by 40 cycles of denaturation at 95 °C for 20 s, annealing at 55 °C for 20 s, and extension at 72 °C for 30 s. Three technical replicates were performed for each sample. Primer sequences are provided in Table 2. GAPDH served as the internal control, and relative gene expression was quantified using the 2^−ΔΔCT method. GraphPad Prism 5 was used for graph creation and statistical analysis, with p-values calculated for each comparison.

Ethics approval and consent to participate

This study was conducted in accordance with the Declaration of Helsinki and approved by the Ethics Committee of The People’s Hospital of Guangxi Zhuang Autonomous Region (2024.03.01, KY-ZC-2024-030).

Consent to participate

Written informed consent was obtained from all participants.

Statistical analysis

All statistical analyses were performed using R software. Differential expression analysis was conducted using the DESeq2 package, and survival analysis was carried out using the survminer package to plot KM survival curves. WGCNA was performed using the WGCNA package, and functional enrichment analysis was conducted with the clusterProfiler package. Cox regression analysis was performed using the R survival package, while LASSO regression analysis was carried out with the glmnet package. Additionally, ROC curves were plotted using the survivalROC package, and nomograms were constructed using the rms package. Immune infiltration cell analysis was conducted using the CIBERSORT algorithm, and the IC50 values for common chemotherapy and molecular targeted drugs were calculated using the pRRophetic package. Correlation analysis was performed using Spearman’s correlation. For pairwise comparisons, the Wilcoxon test was applied, and P < 0.05 was considered statistically significant.

Identification of DEGs and characterization of the correlation between MFRGs and CRC

In TCGA-CRC samples, a total of 3,310 DEGs were identified, comprising 1,669 upregulated and 1,641 downregulated genes (Figs. 1A–1B). To investigate the relationship between MFRGs and CRC, differential analysis revealed that the MFRG score was significantly lower in CRC samples compared to normal samples (Fig. 1C). KM survival analysis showed that groups with high MFRG scores exhibited improved survival probabilities, suggesting that MFRGs may influence the survival of patients with CRC (Fig. 1D). Further analysis of clinical characteristics demonstrated significant differences in MFRG scores across various factors, including race, M stage (M0 vs M1), and N stage (N0 vs N1, N0 vs N2) (Fig. 1E).

Acquisition of 49 DE-MFRGs

Sample clustering and phenotypic trait heatmaps showed no outlier samples (Fig. 2A). A soft-threshold of 7, with an R² value of 0.85, was selected to construct a scale-free network (Fig. 2B). Following average hierarchical clustering and dynamic tree clipping, 16 modules were identified (Fig. 2C). The correlation heatmap revealed that the MEgreen (Cor = 0.67), MEturquoise (Cor = 0.54), MEred (Cor = −0.52), and MEsalmon (Cor = −0.52) modules were significantly associated with CRC (P < 0.05) (Fig. 2D). A total of 917 genes, 3,143 genes, 907 genes, and 302 genes were derived from these modules, resulting in 1,952 hub module genes (GS = 0.4, MM = 0.4) (Figs. 2E–2H). By intersecting the 3,310 DEGs with the 1,952 hub genes, 49 DE-MFRGs were identified (Fig. 2I).

Biological functions and signaling pathways involved in 49 DE-MFRGs

Functional enrichment analysis of the 49 DE-MFRGs revealed 15 significant terms in GO analysis, including caspase binding, helicase activity, kinetochore formation, and cell division site organization (Figs. 3A–3C). KEGG pathway analysis identified three key pathways, such as pentose and glucuronate interconversion and nucleotide sugar biosynthesis (P < 0.05) (Fig. 3D).

Identification of hub gene and risk model construction

The prognostic value of the 49 DE-MFRGs was assessed by constructing a risk model. Initially, univariate Cox analysis was performed on the 49 DE-MFRGs in 606 CRC samples, yielding 10 genes associated with survival (P < 0.05). Of these, HIG1 hypoxia inducible domain family member 1A (HIGD1A), diaphanous related formin 3 (DIAPH3), coiled-coil domain containing 68 (CCDC68) and family with sequence similarity 151 member A (FAM151A) were identified as protective factors, while keratin associated protein 5-1 (KRTAP5.1), heat shock transcription factor 4 (HSF4), melanocortin 1 receptor (MC1R), Zinc finger protein 692 (ZNF692), Lck interacting transmembrane adaptor 1 (LIME1) and tweety family member 3 (TTYH3) were risk factors (Table 3 and Fig. 4A). LASSO analysis, using the minimum lambda value, selected seven genes (HIGD1A, DIAPH3, CCDC68, FAM151A, HSF4, MC1R and TTYH3) to construct the optimal risk model (Figs. 4B–4C). The PH assumption test for these seven genes showed P-values greater than 0.05, and a forest plot from multivariate Cox analysis was generated (Fig. 4D). Finally, stepwise regression analysis identified three hub genes for the construction of a prognostic risk model: CCDC68 (HR = 0.74, 95% CI [0.56–097], P = 0.03) and FAM151A (HR = 0.38 95% CI [0.20–0.71], P = 0.002) were protective factors, and MC1R (HR = 2.22, 95% CI [1.56–3.17], P < 0.001) was a risk factor (Table 4 and Fig. 4E).

Risk model had good predictability for CRC patients

In TCGA-CRC, the high- and low-risk groups divided based on the median of the risk score were calculated, with results showing a significant increase in the number of deaths as the risk scores rose (Figs. 5A–5B). KM analysis revealed that patients in the high-risk group had a lower survival rate (Fig. 5C). ROC analysis showed area under the curve (AUC) values of 0.631, 0.626, and 0.663 at 1-, 2-, and 3-year intervals, respectively (Fig. 5D), indicating that the risk model demonstrated reasonable predictive accuracy. The risk model was then validated in the GSE103479 dataset, with results consistent with TCGA-CRC, confirming that high-risk patients with CRC had worse survival and prognosis (Figs. 5E–5H).

RiskScore, age, and N/M stages were independent prognostic factors for CRC

Further investigation of the correlation between clinicopathologic characteristics and the risk model in TCGA-CRC samples revealed that the risk score, along with age and T/N/M stages, showed significant association with the risk model in univariate Cox analysis (P < 0.05) (Fig. 6A), all of which satisfied the PH assumption (P > 0.05). Multivariate Cox analysis demonstrated that, except for T stage, the remaining clinicopathologic characteristics were significantly correlated with the risk model (P < 0.05) (Fig. 6B). A nomogram was developed to visualize the clinical characteristics, showing that higher scores correlated with increased mortality. The validity of the nomogram was confirmed by calibration curves, with an AUC exceeding 0.6, affirming its strong predictive ability for the survival of patients with CRC (Figs. 6C–6E). Further analysis of the risk scores and clinicopathologic features revealed significant differences in risk scores across T/N/M stages (Fig. 6F).

Enrichment pathways of high- and low-risk groups and their effects in the immune micro environment

GSEA identified the top six signaling pathways significantly associated with the risk model, including cardiac muscle contraction, ascorbate and aldarate metabolism, pentose and glucuronate interconversion, retinol metabolism, extracellular matrix (ECM)-receptor interaction, and dilated cardiomyopathy (Fig. 7A).

Differences in the immune microenvironment between the high- and low-risk groups in CRC were also explored. The proportional distribution of 22 immune cell types is shown in Fig. 7B. Notably, except for M0 macrophages, resting dendritic cells, memory CD4+ T cells, activated memory CD4+ T cells, and plasma cells, the low-risk group exhibited more significant differences in immune cell proportions (Fig. 7C), indicating that the immune system of low-risk group patients was more active and more favorable for controlling tumor progression. Furthermore, 16 immune checkpoints showed significant differences between the high- and low-risk groups, including HHLA2, TNFRSF4, CD160, TNFSF4, and TIMP3 (Fig. 7D). These immune checkpoints may be involved in the activation or inhibition of immune cells, suggesting differences in immune treatment responses between the risk groups. In addition, the TIDE score and immune exclusion score of the high-risk group were significantly higher than those of the low-risk group (Fig. 7E). These two scores are more meaningful in evaluating the effectiveness of immunotherapy. The high-risk group patients, due to their higher scores, were associated with poorer prognosis, including shorter survival and higher recurrence rates. The higher the TIDE score, the poorer the response to immune checkpoint inhibitors and the shorter the survival time.

Fifty-one common chemotherapeutic drugs were effective in high-risk patients for CRC

The IC50 values of 71 chemotherapeutic drugs showed significant differences between the two risk group samples, with 51 drugs demonstrating better responses in the high-risk group patients (IC50 values lower in the high-risk group than in the low-risk group) (Figs. 8–9), while the remaining 20 drugs performed better in the low-risk group (Fig. 10). It was worth noting that 20 of these drugs had already been clinically applied (Table 5). This suggested a potential significant difference in chemotherapy drug sensitivity between the high-risk and low-risk groups. Spearman correlation analysis revealed that FAM151A was positively correlated with AG.014699, camptothecin, NVP.BEZ235, and TW.37 (r > 0.4), while CCDC68 was negatively correlated with BMS.708163 and X681640 (r < −0.4) (Fig. 11). These findings indicate that the expression levels of certain genes may be associated with the sensitivity to specific drugs, suggesting that gene expression levels could influence tumor cell response to these drugs, or the drugs could impact gene-related functions.

Three hub genes linked in other diseases

Pan-cancer analysis showed that CCDC68, FAM151A and MC1R were significantly different in thyroid cancer (THCA), renal clear cell carcinoma (RCC) and lung squamous cell carcinoma (LUSC) etc. (Fig. 12). In TCGA-CRC, the differential expression between CCDC68 and FAM151A were more significant in the normal group, whereas MC1R had opposite results to the above two genes.

Validation of expression of 3 biomarkers

In TCGA-CRC, the differential expression of CCDC68 and FAM151A was more significant in the normal tissue samples, whereas MC1R exhibited opposite expression patterns (Fig. 13A). RT-qPCR validation of hub gene expression in clinical CRC and normal tissues confirmed these findings. CCDC68 and FAM151A were more highly expressed in normal samples compared to CRC samples, while MC1R expression was higher in CRC samples, consistent with the dataset results (Figs. 13B–13D). Notably, the expression of FAM151A (P < 0.0001) and MC1R (P < 0.0001) exhibited significant differences between clinical CRC and normal samples.