RNA sequencing and single nucleotide variation (SNV) data were extracted from The Cancer Genome Atlas (TCGA) database (https://gdc.cancer.gov/). In addition, publicly available data on copy number variation (CNV), methylation, and clinical characteristics of patients were collected. To perform differential expression analysis of MCM family members in various types of cancer, normalized and batch-corrected RSEM mRNA expression data were used. Analysis was conducted on over 10 paired tumours and normal samples from TCGA. Clinical data were obtained from tumour samples of nine cancer types and used to investigate gene expression changes associated with different subtypes. The gene set variation analysis (GSVA) score was used as a measure of gene set expression since it is positively correlated with gene expression. It was calculated using the GSVA R package (version 4.1.1, R Foundation for Statistical Computing, Vienna, Austria). To compare GSVA scores between groups, the Wilcoxon test was used in cases where the number of subtype groups was 2, and ANOVA was used in cases where the number of subtype groups exceeded 2. GSEA was implemented using the fgsea R package (version 4.1.1, R Foundation for Statistical Computing, Vienna, Austria) to determine the overexpression levels of a gene set at the top or bottom of the list of genes. This analysis was based on the gene expression fold change (FC) values between the tumour samples and normal samples.

The GSCALite platform was used for the genome-wide analysis of the MCM genes [17]. Differential expression of MCM family members was evaluated using RSEM-normalised RNA-seq data. A heatmap and waterfall plot were constructed to visualise the SNV and CNV frequencies of the MCM genes. The CNV data were obtained from 11,495 samples in the TCGA database, while GISTIC2.0 was used to identify regions of significant amplification or deletion in various patient groups [18].

The DNA methylation data (Illumina Human Methylation 450k level 3) were extracted from the TCGA database for 14 cancer types. More than 10 pairs of tumours and adjacent normal tissues were analysed. Multiple methylation sites were typically present in a single gene, each with its tag storing the methylation level.

In this study, we collected clinical data from 33 different types of cancers for the purpose of conducting survival analysis. To integrate methylation and clinical survival data, sample barcodes were used, and tumour samples were categorized into high- and low-methylation groups based on the median methylation level. In this study, we utilized the survival package in R (version 4.1.1, R Foundation for Statistical Computing, Vienna, Austria) to build Cox proportional hazards models and conduct log-rank tests on MCM genes across various types of cancers. A p-value of less than 0.05 was deemed statistically significant.

Pathway activity scores (PAS) for cancer-related pathways were calculated using RPPA data extracted from the TCPA (The Cancer Proteome Atlas, http://www.tcpaportal.org) database. In order to estimate differences in PASs between groups, the Student’s t-test was utilized, with p-values being adjusted through the False Discovery Rate (FDR) method. Any FDRs that were equal to or less than 0.05 were considered to indicate significant differences [19, 20].

The IC${}_{50}$ values of the small molecules and corresponding mRNA expression data were sourced from the Genomics of Cancer Drug Sensitivity (GDSC) and Genome Therapy Response Portal (CTRP). The drugs were ranked according to the correlation coefficients of the screened genes and the FDR comprehensive levels. Drug–gene pairs with absolute correlation coefficients of $>$0.1 and FDRs of 0.05 were retained, and each pair was assigned a score by multiplying their –log${}_{10}$FDRs and absolute correlation coefficients.

The study utilized the immune infiltration and GSVA scoring modules to investigate the correlation between gene expression and immune cell infiltration. Gene set expression was assessed by estimating GSVA scores, which were positively correlated with gene set expression.

The study utilized Student’s t-test to determine variations in the MCM family expression and methylation levels in tumour tissues compared to their levels in corresponding normal tissues. The log-rank test was utilized to compare survival curves, while Pearson or Spearman analyses were utilized to estimate correlation coefficients. A p-value of less than 0.05 was deemed statistically significant. All statistical analyses were conducted using SangerBox database (http://www.sangerbox.com/, Sanger Box 3.0; Hangzhou, China) [21] and ChiPlot (https://www.chiplot.online/, Shantou, China), free online platforms for data analysis (*, p 0.05; **, p 0.01; ***, p 0.001; ****, p 0.0001).

Paired tumour and normal tissue gene expression data from 33 types of cancers in the TCGA database were used to analyse the differential expression of the MCM genes. The results revealed that MCM2, MCM3, MCM4, MCM5, MCM6, and MCM7 were upregulated in most tumour tissues, whereas MCM8 and MCM9 were downregulated in certain types of tumours. Specifically, MCM8 was downregulated in kidney renal clear cell carcinoma (KIRC) and thyroid carcinoma (THCA), whereas MCM9 was downregulated in kidney chromophobe (KICH), KIRC, and kidney renal papillary cell carcinomas (KIRP) (Fig. 1A). GSVA was used to determine the gene set expression. The GSVA scores for bladder urothelial carcinoma (BLCA), breast invasive carcinoma (BRCA), colon adenocarcinoma (COAD), oesophageal carcinoma (ESCA), head and neck squamous cell carcinoma (HNSC), KICH, KIRC, KIRP, liver hepatocellular carcinoma (LIHC), lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC), stomach adenocarcinoma (STAD), and THCA were significantly higher than in the corresponding adjacent normal tissues (Fig. 1B). The expressions of MCM2–8 have been found to be significantly correlated with BRCA subtypes. Notably, there was a significant difference in the expressions of MCM2–8 proteins between the subtypes of luminal A and luminal B, thereby suggesting that they could potentially be used as biomarkers to differentiate between the two subtypes (Fig. 1C, Supplementary Fig. 1). An assessment of the standardised enrichment scores (NESs) revealed that the MCM gene set was significantly enriched in several types of cancers, including HNSC, ESCA, COAD, LUSC, LIHC, BLCA, STAD, BRCA, and LUAD (Fig. 1D), especially in HNSC and ESCA (Fig. 1E,F). These findings indicate that irregular expression of the MCM genes plays a role in the onset of multiple forms of cancer.

Fig. 1.

Comparison of the Minichromosome Maintenance Proteins (MCM) family mRNA expression in various types of cancer. (A) Differential expression of MCM genes between 14 paired normal and tumour tissues. Purple to red represents the fold change between tumour tissues and normal tissues. Red dots indicate that gene expression is higher in tumour tissues than in normal tissues, and blue dots indicate that gene expression is lower in tumour tissues than in normal tissues. The size of the dots indicates significance based on the False Discovery Rate (FDR), while the FDR value of each dot is shown on the right. (B) The gene set variation analysis (GSVA) score was calculated by comparing the molecular characterization maps (MCMs) between 14 paired normal and tumour tissues. (C) Relationship between the MCM gene expressions and subtypes. FDR values are represented by the bubble’s colours and sizes in the graph. Rows represent gene set symbols, whereas columns represent selected cancer types. Significance based on FDR values is indicated by the colour of the bubbles, with white indicating low significance and red indicating high significance. Bubble size positively correlates with significance. (D) Summary of the enrichment scores (ES) of MCM genes in the selected cancers. The score indicates the extent to which a set of genes is overrepresented at either the top (ES $>$0) or bottom (ES 0) of a gene ranking list. (E,F) MCM gene sets were found to be significantly enriched in head and neck squamous cell carcinoma (HNSC) and oesophageal carcinoma (ESCA). In the plot, the running ES for the gene set is displayed as the analysis descends the ranked list. DEGs, Differentially expressed genes; KICH, kidney chromophobe; PRAD, prostate adenocarcinoma.

To identify the various genetic alterations in the MCM genes, the frequency and types of SNVs were evaluated across the 33 types of cancers. Uterine corpus endometrial carcinoma (UCEC) had the highest SNV frequency (198%), followed by skin cutaneous melanoma (SKCM) (110%), STAD (77%), LUAD (67%), COAD (66%), BLCA (63%), LUSC (56%), HNSC (40%), BRCA (32%), glioblastoma multiforme (GBM) (26%), LIHC (25%), cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC) (22%), rectum adenocarcinoma (READ) (19%), OV (15%), KIRC (12%), prostate adenocarcinoma (PRAD) (12%), and brain low-grade glioma (LGG) (12%). The SNV frequency was 10% among the remaining 16 cancer types, and notably, no mutations were observed in the MCM genes in thymoma (THYM) (Fig. 2A). Conversely, genes with high SNV frequency included MCM4 (21%), MCM6 (19%), MCM7 (18%), MCM2 (17%), MCM3 (16%), MCM8 (16%), MCM5 (15%), and MCM9 (12%) (Fig. 2B). Furthermore, we evaluated the various mutation frequencies in the MCM genes. The study found that missense mutations were the most frequently identified mutation, with MCM4 and MCM6 having the highest mutation frequency. In particular, missense mutations were the most common type observed. Thereafter, a pie chart was created to provide a visual representation of the relationship between the different CNV types and mRNA expressions of the MCM members. Heterozygous amplification and heterozygous deletion were found to be the primary CNV types (Fig. 2C). Multiple cancer types exhibited heterozygous amplification of MCM2, MCM3, MCM4, MCM5, and MCM6. Specifically, a heterozygous amplification frequency of $>$25% was observed for MCM2 in BLCA, CESC, ESCA, HNSC, KIRP, LUSC, OV, STAD, and uterine carcinosarcoma (UCS); MCM3 in ESCA, LIHC, LUAD, LUSC, OV, READ, SKCM, UCS and uveal melanoma (UVM); MCM4 in adrenocortical carcinoma (ACC), BLCA, BRCA, CESC, COAD, ESCA, HNSC, KICH, LIHC, LUAD, LUSC, OV, READ, sarcoma (SARC), SKCM, STAD, testicular germ cell tumours (TGCT), UCEC, UCS and UVM; MCM5 in KICH, LUSC, SKCM, and TGCT; MCM6 in LUAD, LUSC, OV, TGCT, and UCS; MCM7 in ACC, BLCA, BRCA, COAD, lymphoid neoplasm diffuse large B-cell lymphoma (DLBC), ESCA, GBM, HNSC, KICH, KIRC, KIRP, LGG, LIHC, LUAD, LUSC, mesothelioma (MESO), OV, pancreatic adenocarcinoma (PAAD), READ, SARC, SKCM, STAD, TGCT, and UCS; MCM8 in ACC, BLCA, BRCA, CESC, cholangiocarcinoma (CHOL), COAD, ESCA, GBM, HNSC, KICH, KIRP, LIHC, LUAD, LUSC, OV, READ, SKCM, STAD, and UCS; MCM9 in UCS. Our research found a notable correlation between the mRNA expression of the MCM members and CNV frequency in the majority of cancer types (p 0.05, Fig. 2D). These findings indicate that the MCM genes are subject to amplification and loss of heterozygosity mutations, which lead to abnormal expression, and ultimately, contribute to the development of cancer.

Fig. 2.

Relationship between the MCM gene expressions and genomic alterations. (A) Heatmap demonstrating the SNV frequency in MCM2–9 genes in tumours. (B) Waterfall plot demonstrating SNV distribution and its classification in MCM genes. (C) Pie chart demonstrating the proportion of different types of copy number variations (CNVs) in each gene across various cancer types. (D) Bubble plots demonstrating the relationship between mRNA expression and CNV frequency. The colour of the bubbles indicates the strength of correlations, with blue bubbles representing negative correlations and red bubbles representing positive correlations. The bubble sizes are proportional to the significance, based on the FDR values. Black borders indicate an FDR of $\le$0.05.

The study found a strong negative correlation between the methylation levels and mRNA expression in most MCM family members. However, this correlation was not observed in MCM7 in CHOL and MCM3 in CESC, DLBC, KICH, KIRC, acute myeloid leukaemia (LAML), OV, PAAD, pheochromocytoma and paraganglioma (PCPG), PRAD, THYM, and UCS (Fig. 3A). Analysis of the differential methylation revealed that the methylation levels of the MCM genes were significantly lower in tumour tissues compared to the matched normal tissues: MCM2 in BLCA, BRCA, ESCA, HNSC, KIRC, KIRP, LIHC, LUAD, LUSC, PAAD, THCA, and UCEC; MCM3 in BLCA, BRCA, LIHC, LUSC, and PRAD; MCM4 in BRCA, ESCA, KIRC, LUAD, and THCA; MCM5 in BLCA, BRCA, COAD, KIRC, LIHC, LUAD, LUSC, PAAD, PRAD, and UCEC; MCM6 in BLCA, BRCA, HNSC, KIRC, LIHC, LUAD, LUSC, THCA, and UCEC; MCM7 in KIRP, LIHC, LUAD, PRAD, and THCA; MCM8 in BRCA and PRAD; MCM9 in BLCA, KIRP, LIHC, LUAD, LUSC, PRAD, and UCEC. However, the methylation levels of the following genes were higher in the tumour tissues than in the matched normal tissues: MCM3 in KIRC, KIRP, and THCA; MCM4 in KIRP; MCM5 in ESCA; MCM6 in PRAD; MCM7 in BRCA and LUSC; MCM8 in KIRC, LUSC, and THCA; MCM9 in BRCA, KIRC, and PAAD (Fig. 3B). Hypomethylation of MCM2–9 in various types of cancers was associated with an increased risk of mortality (Fig. 3C). The study found that hypermethylation of MCM2, MCM3, MCM7, and MCM8 was a significant risk factor for survival in certain types of cancers, including UVM, LGG, CESC, and ACC (Fig. 3D–H). Altogether, these results suggest that abnormal DNA methylation regulates abnormal expression of MCMs, thereby influencing tumour progression.

Fig. 3.

Methylation alterations in MCM2–9 genes and overall survival (OS). (A) Correlation between the mRNA expression and methylation levels of MCM genes. Negative correlations are represented by blue bubbles, whereas positive correlations are represented by red bubbles. The intensity of the colour reflects the strength of the correlation, with darker colours indicating a stronger correlation. The size of the bubbles is positively correlated with significance, based on FDR values. Black borders indicate an FDR of $\le$0.05. (B) Differential methylation of MCM genes between 14 paired normal and tumour tissues. The colour and size of the bubbles represent fold change values and FDRs, respectively. Rows represent gene symbols, whereas columns represent the selected cancer types. (C) Differences in overall survival between hypermethylated and hypomethylated groups of MCM genes in specific cancers. The plot displays hazard ratios and p-values (Cox regression analysis) for various cancer types and gene symbols. The colour of the bubbles indicates the hazard ratio from low to high, whereas the size of the bubbles represents a positive correlation with the p-value. Bubbles with black borders signify p-values of $\le$0.05. (D–H) KM plots (version 4.1.1, R Foundation for Statistical Computing, Vienna, Austria) demonstrating differences in overall survival between patients with low and high hypermethylation groups for MCM2 in uveal melanoma (UVM), MCM3 in brain low-grade glioma (LGG), MCM7 in cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), and MCM8 in adrenocortical carcinoma (ACC) and CESC.

The survival analysis results indicate that in most cancer types, a high expression of MCM2–9 was linked to a poor prognosis. However, a high expression of certain MCM genes was associated with a better prognosis in specific cancer types: MCM2 in CESC, STAD, and UVM; MCM3 in CESC, OV, STAD, and THCA; MCM4 in READ and THYM; MCM5 in CESC and THCA; MCM6 in CESC, READ, and THYM; MCM7 in DLBC and THCA; MCM9 in THYM (Fig. 4A–C). Various tumours were found to be influenced by the MCM2–9 genes, in terms of their pathological staging. Specifically, MCM4, MCM6, and MCM2 were associated with KIRP; MCM7, MCM5, MCM6, and MCM2 were associated with KIRC; MCM3 and MCM2 were associated with BRCA; MCM5 and MCM6 were associated with THCA; MCM7 was associated with KICH; MCM8 was associated with SKCM; MCM3 was associated with TGCT (Fig. 4D). According to the findings for KICH, a positive correlation existed between MCM2 expression and the pathological stage. This suggests that tumours with higher stages tend to exhibit higher levels of MCM2 expression (Fig. 4E). Furthermore, the MCM2–9 genes exhibited varying levels of expression in different clinical stages of tumours. Specifically, MCM4, MCM6, and MCM2 were differentially expressed in different clinical stages of KIRP (Fig. 4F). While the expression of MCM3 increased with the clinical stage in TGCT (Fig. 4G). This study suggests that the abnormal expression of MCM genes is closely linked to cancer prognosis.

Fig. 4.

Relationship between MCM gene expressions and overall survival (OS) in pan-cancer. (A) Heatmap demonstrating the univariate Cox proportional-hazards regression analysis –log${}_{10}$ value (p-value). The colour of the bubbles from blue to red represents the low to high hazard ratios, respectively, whereas the size of the bubbles positively correlates to the p-value. Bubbles with black borders have p-values of $\le$0.05. (B,C) KM plots demonstrate differences in overall survival between the low and high MCM expression groups in UVM and OV. Statistical analysis was performed using the log-rank test, with p-values of 0.05 indicating significance. (D) Differences in the mRNA expression of MCM genes among different pathological stages of specific cancers. (E) Boxplot demonstrating MCM2 expression in different pathological stages of KICH. (F) Differences in the mRNA expression of MCM genes among different clinical stages in specific cancers. (G) Boxplot demonstrating MCM3 expression in different clinical stages of TGCT. *, p 0.05; **, p 0.01; ***, p 0.001; NS, no significance.

This study investigated the correlation between MCM gene expressions and immune cell infiltration in the tumour microenvironment (TME). MCM2 expression significantly correlated with the immune, stromal, and microenvironmental scores. A negative correlation was observed between the immune scores and MCM2 expression in cancers, including GBM, UCEC, ESCA, stomach and oesophageal carcinoma (STES), SARC, KIRP, LUSC, high-risk Wilms tumour (WT), neuroblastoma (NB), TGCT, and PCPG (Fig. 5A). In this study, we investigated the relationship between the infiltration levels of 24 immune cell types and GSVA scores. Our findings suggest the existence of two distinct clusters of cancer types based on this correlation (Fig. 5B). The two clusters exhibited different patterns of immune cell infiltration. The MCM gene expressions positively correlated with the infiltration levels of the immunosuppressive cells in cancers, including the natural regulatory T (nTreg) cells, induced regulatory T (iTreg) cells, exhausted, dendritic cells (DCs), and macrophages. However, the expression of the MCMs was significantly negatively correlated with the infiltration levels of the immune effector cells, including the natural killer (NK) cells, CD8${}^{+}$ T cells, follicular helper T (Tfh) cells and CD4${}^{+}$ T cells, in cancers, such as in TGCT, THCA, HNSC, LUAD, STAD, CESC, PCPG, THYM, GBM, and KIRP. These findings suggest that the upregulation of the MCM gene expressions is associated with the immune microenvironment of these tumours. Therefore, targeting the MCM family could have significant implications for improving the effectiveness of immunotherapy in cancer patients.

Fig. 5.

Relationship between the MCM gene expressions and the tumour microenvironment. (A) Heatmap demonstrating the correlation between MCM2 expression and immune, stromal, and microenvironmental scores. (B) Relationship between the MCM gene expressions and immune cell infiltration. The heatmap summarises the significant correlations based on p-values and FDRs from Spearman correlation analysis of GSVA scores and immune cell infiltration. Blue represents the negative correlations, whereas red represents positive correlations. Note: *, p 0.05; **, p 0.01; ***, p 0.001; ****, p 0.0001; #, FDR 0.05.

To investigate the impact of the MCM genes on drug efficacy, the study analysed the correlation between gene expression and drug sensitivity using data obtained from the GDSC and CTRP databases. The results indicated that the increased expression of the MCM2–9 genes was linked to heightened resistance to drugs, such as 17-AAG, RDEA119, trametinib, and selumetinib (Fig. 6A). Moreover, the study found that there was a negative correlation between the expressions of MCM2–9 and drug sensitivity. This correlation was validated by analysing the IC${}_{50}$ values of the drugs (Fig. 6B). These findings suggest that the overexpression of the MCMs results in a resistance to chemotherapy and targeted drug therapy, meaning that MCMs are useful predictors of drug efficacy and potential therapeutic targets for cancer.

Fig. 6.

Relationship between the MCM gene expressions and drug sensitivity. (A) Correlation between the expressions of the MCM genes and drug sensitivity based on Genome Therapy Response Portal (CTRP) data. (B) Correlation between the expressions of the MCM genes and drug sensitivity based on GDSC data. In the scatterplot, blue bubbles indicate negative correlations and red bubbles indicate positive correlations, while the bubble size positively correlates to significance based on FDR values. Darker colours indicate stronger correlations. Additionally, black borders indicate FDR values of $\le$0.05.

The miRNA–gene network analysis revealed that all eight MCM genes were regulated by multiple miRNAs (Fig. 7A). The expressions of the MCM genes were regulated by a complex network of miRNAs, which play an essential role in tumour development and progression. Pathway enrichment analysis revealed that the MCM2–9 genes were associated with several pathways, including the TSC/mTOR signalling, cell cycle-related, PI3K/AKT signalling, RTK signalling, RAS/MAPK signalling, and apoptosis-related pathways (Fig. 7B). The MCM genes had a positive impact on specific pathways in several cancers, whereby MCM7 and MCM5 were primarily involved in the activation of Androgen Receptor (AR) and Epithelial–Mesenchymal Transition (EMT), respectively (Fig. 7C,D). The GSVA demonstrated that the expressions of the MCM gene sets were linked to cell cycle regulation and apoptosis in cancer (Fig. 7E). This finding suggests that the MCMs have a significant role in controlling cancer-related pathways.

Fig. 7.

Potential molecular mechanisms underlying altered MCM gene expressions. (A) The miRNA regulatory network is representative of miRNAs and target genes, whereby nodes represent miRNAs or target genes, edges represent the miRNA-to-gene conversion regulation, and edge widths indicate the absolute correlation coefficients. (B) The combined percentage of the effects of the MCM complex on pathway activity. (C) Pie chart demonstrating the proportion of genes contributing to the pathway activity across 32 cancer types. (D) A network demonstrating the connection between genes and pathways using straight lines. Activation is represented by solid lines, whereas inhibition is denoted by dashed lines. The line colours are used to indicate different cancer types. (E) Association between GSVA scores and cancer-related pathway activity. *, p-value $\le$ 0.05; #, FDR $\le$0.05.