Copper Metabolism-Related Genes as Biomarkers in Colon Adenoma and Cancer

Introduction

Colorectal cancer (CRC) ranks as the third most frequently diagnosed malignant tumor globally and remains a major contributor to cancer-related deaths.¹ As of 2019, the global incidence of CRC has risen to 2.19 million cases, with Chinese patients comprising approximately 30% of the worldwide CRC population.² According to recent projections, the annual number of new CRC cases is expected to rise to 2.5 million by 2035, posing a growing public health challenge worldwide.³ The treatment of CRC imposes a significant financial burden, with direct medical costs estimated at $106.4 billion, accounting for 62.4% of the total economic impact.⁴ The rapid increase in cases and high economic costs highlight the urgent need for effective prevention, early detection, and treatment. The progression from normal colonic epithelium to adenoma and invasive carcinoma, driven by acquired molecular events, usually takes 5–10 years.^5,6 Colonic adenoma (CA) is a benign epithelial tumor arising in the colon’s inner lining, marked by abnormal cell proliferation and the formation of small masses. Approximately 10% of CA cases may progress to CRC, influenced by the accumulation of gene mutations, epigenetic abnormalities, and environmental factors.^7,8 Colon cancer progression is closely linked to various genetic mutations, with mutations in the APC, KRAS, and TP53 genes being particularly prevalent in CRC.⁹ APC mutations lead to the over-activation of the Wnt signaling pathway, a critical driver of tumorigenesis.¹⁰ Similarly, mutations in the KRAS gene facilitate uncontrolled tumor cell proliferation.¹¹ Furthermore, TP53 mutations compromise the cell’s ability to respond to DNA damage, thereby accelerating cancer cell development and progression.¹² Mutations in DNA repair genes (eg, MLH1, MSH2, MSH6) are strongly linked to hereditary nonpolyposis colorectal cancer (HNPCC). They disrupt DNA repair, increasing mutation accumulation and driving CRC development.¹³ Epigenetic alterations also play a pivotal role in CRC progression. For instance, mutations and reduced expression of the histone H3K36 trimethyltransferase SETD2 are linked to the advancement of CRC.¹⁴ Additionally, aberrant DNA methylation, particularly the hypermethylation of promoters associated with microsatellite instability (MSI), leads to reduced or absent expression of DNA mismatch repair proteins, further contributing to CRC pathogenesis.¹⁵ Meanwhile, environmental factors cannot be ignored. Unhealthy dietary habits, obesity, and long - term social stress can cause DNA damage. Once the DNA was damaged and the repair mechanism is compromised due to epigenetic changes, the risk of developing CRC significantly increases.¹⁶ Evidence suggests that the detection and excision of adenomas substantially decrease CRC incidence and mortality.¹⁷ Therefore, early CA detection is critical. Although colonoscopy remains effective for identifying and excising adenomas before their progression to CRC, it has notable limitations, including a missed diagnosis rate of up to 25%.¹⁸ Stool-based screening, commonly employed for CRC detection, lacks sensitivity for adenomas, making it unsuitable for early detection of these lesions.^8,19 This underscores the urgent need to identify novel characteristic genes associated with CA.

Copper (Cu) is a vital trace element required for the growth and development of all eukaryotes.²⁰ It plays a critical role in redox enzymatic activities that are essential for key metabolic processes, signaling pathways, and various biological functions.²¹ The link between disrupted Cu metabolism and the pathogenesis of various malignancies has attracted significant research attention. Studies indicate that inflammatory cytokines, such as interleukin (IL)-17, enhance Cu uptake through inflammatory responses, thereby fostering colon tumor development.²² Furthermore, evidence shows that Cu levels in colon tumor tissues are significantly elevated compared to normal colon mucosa.²³ Studies have shown that copper metabolism-related genes (CMRGs) are associated with immune infiltration in colon adenocarcinoma (COAD). They can be utilized to classify COAD subtypes and distinguish patients with different molecular characteristics, serving as valuable biomarkers for immunotherapy in COAD.²⁴ This highlights Cu metabolism’s critical role in colon tumor initiation and progression. However, research on CMRGs in the context of colon cancer development remains limited.

Bioinformatics analysis serves as a powerful tool for predicting molecular pathways and gene interactions, and it has been extensively used to identify novel candidate genes and pathways across various cancer types.²⁵ By leveraging computational methods to analyze extensive genomic and transcriptomic data, bioinformatics has facilitated the identification of novel candidate genes and pathways associated with the development of colon carcinogenesis.^26,27 This study aims to conduct a comprehensive bioinformatics analysis to identify signature genes that can effectively distinguish between colon adenomas (CAs), colorectal cancers (CRCs), and normal colons (NCs). Based on these identified genes, we expect to develop a robust classification model that can accurately differentiate between CA, CRC, and NC, thereby enhancing both the specificity and sensitivity of colorectal disease diagnosis.Furthermore, the study will incorporate single-cell data analysis from NC, CA, and CRC patient samples, with the aim of elucidating the expression patterns of these signature genes across various cell types, including fibroblasts, lymphocytes, cuprocytes, B-cells, and macrophages. This analysis is expected to provide valuable insights into the cellular mechanisms that underpin disease progression and to shed light on the roles of these genes within the tumor microenvironment. By identifying these key genes and their interactions, the study seeks to offer novel perspectives on the molecular mechanisms driving CRC progression, thereby advancing the field of colorectal cancer research and informing future diagnostic strategies.

Materials and Methods

Data Source

The RNA sequencing (RNA-seq) datasets GSE117606, GSE20916 and GSE41657, along with the single-cell RNA-seq (scRNA-seq) dataset GSE201348, were retrieved from the Gene Expression Omnibus (GEO) database (available at https://www.ncbi.nlm.nih.gov/geo/). The GSE117606 dataset, detected using the GPL25373 platform (Affymetrix HT HG-U133+ PM Array Plate), includes tissue samples from 65 normal controls (NC), 69 patients with CA, and 74 with CRC. The GSE20916 dataset, detected on the GPL570 platform (Affymetrix Human Genome U133 Plus 2.0 Array), consists of tissue samples from 44 NCs, 55 patients with CA, and 46 patients with CRC. The GSE41657 dataset, detected using the GPL6480 platform (Agilent-014850 Whole Human Genome Microarray 4x44K G4112F (Probe Name version)), includes tissue samples from 12 NC, 51 patients with CA, and 25 with CRC. The GSE201348 dataset, processed with the GPL24676 platform (Illumina NovaSeq 6000), contains tissue samples from 8 NCs, 43 patients with CA, five patients with CRC, and 16 unaffected CRC individuals. A total of 1,874 CMRGs were sourced from the GeneCards database (accessible at https://www.genecards.org/) using the search term “Cu metabolism” under the criteria: category = protein-coding and relevance score ≥1. Furthermore, the Cancer Genome Atlas (TCGA)-CRC datasets, comprising 598 CRC patient tissue samples, were accessed from the TCGA database (available at https://gdc.cancer.gov/) for feature gene survival analysis.

Differential Expression Analysis

Differential expression analysis was the process of identified genes with significant differences in expression between sample groups. To obtained the genes with significant differences in the samples, the differentially expressed genes (DEGs) between CA and NC samples in the GSE117606 dataset were identified using the limma package, applying thresholds of |logFC| > 0.5 and adj. p < 0.05.²⁸ The top 50 DEGs1 were visualized using the pheatmap package (version 1.0.12).²⁹ Differentially expressed CMRGs (DE-CMRGs1) were subsequently identified by intersecting DEGs1 with CMRGs through the VennDiagram package (version 1.6.2).³⁰ Similarly, DEGs2 between CRC and CA samples in the GSE117606 dataset were also determined using the same limma parameters, and DE-CMRGs2 were obtained by overlapping DEGs2 with CMRGs. To explore common functions and relevant pathways, ingenuity pathway analysis (IPA) software was employed to conduct classical pathway enrichment analysis for DE-CMRGs1 and DE-CMRGs2 separately.

Diagnostic Model Construction and Validation

In order to comprehensively validated the diagnostic performance and prediction ability of the model and biomarkers, multiple algorithms were integrated for key marker screening. Given the strong reliability of integrating multiple algorithms for the identification of key markers,^31,32 five machine learning techniques—random forest (RF), support vector machines (SVM), logistic regression (LR), gradient boosting decision tree (GBDT), and artificial neural network (ANN)—were applied to the overlapping genes of DE-CMRGs1 and DE-CMRGs2 to develop diagnostic models for distinguishing between samples from patients with CA, CRC, and NC. The predictive performance of each model, in terms of the area under the curve (AUC) value, was assessed using the pROC package (version 1.17.0.1).³³ The model yielding the highest AUC was selected as the optimal model, and the top 10 most important variables within this model were identified and ranked. Subsequently, biomarkers were refined by overlapping the top 10 variables from the three best-performing models, and the diagnostic value of each individual feature gene was analyzed across different subgroups (CA vs NC; CRC vs CA; CRC vs NC). Additionally, the collective diagnostic accuracy of all biomarkers was evaluated using a logistic regression model through the rms package. The expression levels of biomarkers among samples from NC, CA, and CRC in the GSE117606 dataset were also examined, and their associations with CRC survival in the TCGA-CRC cohort were investigated using Kaplan-Meier (K-M) analysis. Optimal expression thresholds for different groups were calculated using the survminer (version 0.4.9) and survival (version 3.3–1) packages. To further facilitate clinical decision-making regarding CA outcomes, a nomogram was constructed using logistic regression, wherein the score for each feature gene was computed and aggregated into a total score to predict the risk of CA. The diagnostic accuracy of the nomogram was evaluated through calibration curves, decision curve analysis (DCA), and clinical impact curves. The expression levels and diagnostic relevance of the biomarkers were further validated in an external dataset (GSE20916 and GSE41657), with the diagnostic performance of individual biomarkers and the full set analyzed in a manner consistent with the procedures used for the GSE117606 and GSE41657dataset.

Single Gene Set Enrichment Analysis (GSEA)

To comprehensively investigate the pathways influenced by the biomarkers and their associated functions, GSEA was performed using the ClusterProfiler package (version 3.18.0).³⁴ Initially, correlation coefficients were calculated through correlation analysis between the biomarkers and other genes. GSEA was then conducted based on the ranking of these correlation coefficients. The results were visualized using GSEA spine plots generated by the Enrichplot package (version 1.10.2).³⁵

Immune Infiltration Analysis

In order to further explored the pathways targeted by biomarkers and the functions, a single-sample GSEA (ssGSEA) was also performed to analyze immune infiltration, specifically assessing the tumor microenvironment (TME) cell composition in samples from NC and patients with CA and CRC. The proportions of immune cells among these three groups were compared using the rank sum test. Additionally, Spearman correlation analysis was employed to explore the relationships between biomarkers and immune cells.

scRNA-Seq Analysis

The main objective of single-cell RNA sequencing (scRNA-seq) analysis was to study the heterogeneity of gene expression at the single-cell level, to reveal cell-to-cell variability, as well as to explore cell development, differentiation, and function.³⁶ To explore the characteristic gene expression patterns during cell state transitions, as well as cellular changes in disease states, we performed scRNA-seq analysis. For the scRNA-seq data, processing began with the CreateSeuratObject function of the Seurat package (version 4.0.5). The data were filtered using minimum criteria of 100 cells and 100 features.³⁷ The percentages of erythrocyte and mitochondrial genes were calculated using the PercentageFeatureSet function. Data normalization was performed using the NormalizeData function, and the FindVariableFeatures function was utilized to identify genes with high intercellular expression variation. To ensure uniformity, the data were further processed using ScaleData, followed by principal component analysis (PCA). The JackStraw and ScoreJackStraw functions were used for linear dimensionality reduction to assess data quality. The information content of each principal component was ranked, and a PCA ElbowPlot was generated. The JackStrawPlot was employed to compare the distribution of p-values for each principal component. Unsupervised cluster analysis was subsequently performed using the Seurat package, with cell types identified through auxiliary annotation using the SingleR function on the Cellmarker website. The expression levels of biomarkers across different categories were calculated, and their relationships with cell categories were visualized using heatmaps and scatter plots.

Construction of the ceRNA Network

The miRWalk database was used to predict miRNAs targeting the biomarkers, with criteria including an energy score below −15 and a binding probability above 0.95. The Starbase database (accessible at https://starbase.sysu.edu.cn/index.php) was employed to predict the targeting relationships between lncRNAs and miRNAs, using parameters such as a clipExpNum greater than 1, a degraExpNum greater than 1, and a gene type of lncRNA. This information was used to construct the ceRNA regulatory network in Cytoscape (version 3.8.2), analyzing the relationships between lncRNAs, miRNAs, and mRNAs.

Potential Therapeutic Drug Prediction

The Drug-Gene Interaction database (DGIdb, accessible at https://dgidb.genome.wustl.edu/) was utilized to predict drugs related to the biomarkers, aiming to identify small-molecule compounds with potential therapeutic effects for patients with CA or CRC. The gene-drug interactions were visualized using Cytoscape, providing a comprehensive overview of potential therapeutic targets.

Expression Verifications of Biomarkers

To further characterize the mRNA and protein expression of the biomarkers in clinical samples, human tissue samples were obtained, including 10 pairs of frozen samples from patients with CA and NC,³⁸ from the Second Affiliated Hospital of Guizhou University of Traditional Chinese Medicine for qPCR experiments. This study was approved by the Ethics Committee of the Second Affiliated Hospital of Guizhou University of Traditional Chinese Medicine, and written informed consent was obtained from all participants. Total RNA was extracted from the 10 pairs of frozen samples, each samples 50mg, and RNA concentration was measured using the NanoPhotometer N50 (Table S1). Reverse transcription was performed with the SureScript First-Strand cDNA Synthesis Kit from Saville. The reverse transcription system was shown in Table S2. The qPCR reactions were carried out in a CFX96 real-time quantitative fluorescence PCR instrument using a system comprising 3 µL of cDNA, 5 µL of 2x Universal Blue SYBR Green qPCR Master Mix, and 1 µL each of upstream and downstream primers. The amplification protocol included an initial denaturation at 95°C for 1 minute, followed by 40 cycles of denaturation at 95°C for 20 seconds, annealing at 55°C for 20 seconds, and extension at 72°C for 30 seconds. Relative gene expression levels were calculated using the 2-ΔΔCt method, with glyceraldehyde 3-phosphate dehydrogenase (GAPDH) as the internal reference gene.The primer sequences for PCR were shown in Table 1. Additionally, the protein levels of the biomarkers in tissues from patients with CRC were assessed using immunohistochemistry (IHC) data available from the Human Protein Atlas (HPA) database (accessible at https://www.proteinatlas.org/).

Table 1 Sequence of Primers

Statistical Analysis

Statistical analysis was performed using R software. Differences in gene expression and immune cell proportions between different groups were analyzed via the Wilcoxon test, while the qPCR expression data was analyzed using a t-test. Statistical significance was set at adj. p < 0.05.

Results

Identification of DE-CMRGs1 and DE-CMRGs2

The differentially expressed genes (DEGs) between samples from patients with CA and NC, a total of 1,343 DEGs1 were identified between samples from patients with CA and NC, with 467 genes up regulated and 876 genes down regulated (Figure 1A). To visualize these findings, a heat map depicting the top 50 DEGs1 was generated (Figure S1A). By intersecting the set of 1,343 DEGs1 with the 1,874 CMRGs, 218 DE-CMRGs were identified, as illustrated in the Venn diagram (Figure 1B). Similarly, between samples from patients with CRC and CA, 982 DEGs were identified, including 641 up regulated genes and 341 down regulated genes (Figures 1C and S1B). By overlapping these 982 DEGs with the 1,874 CMRGs, 185 DE-CMRGs2 were identified, as shown in the Venn diagram (Figure 1D).

Figure 1 Differentially expressed copper metabolism-related genes (DE-CMRGs) were identified by overlapping the differentially expressed genes (DEGs1) between colon adenoma (CA) and normal control (NC) with copper metabolism-related genes (CMRGs), as well as overlapping the DEGs2 between CRC and CA with CMRGs. The analysis yielded the following results: (A) A volcano plot displaying 1,343 DEGs1 between samples from patients with CA and NC (|logFC| > 0.5 and p < 0.05), red for up-regulated genes, blue for down-regulated genes; (B) Venn diagram illustrating the 218 DE-CMRGs1; (C) A volcano plot showing 982 DEGs2 between samples from patients with CRC and CA (|logFC| > 0.5 and p < 0.05), red for up-regulated genes, blue for down-regulated genes; (D) Venn diagram depicting the 185 DE-CMRGs2.

Identification of DE-CMRGs1 and DE-CMRGs2 Enrichment Pathways

IPA was performed to investigate the functions and associated pathways of DE-CMRGs1 and DE-CMRGs2. The analysis revealed that DE-CMRGs1 were involved in the activation of aryl hydrocarbon receptor signaling while inhibiting the TME pathway, oncostatin M signaling, and LXR/RXR activation. In contrast, DE-CMRGs2 were found to activate the TME pathway, oncostatin M signaling, pulmonary fibrosis idiopathic signaling, the role of chondrocytes in the rheumatoid arthritis signaling pathway, and hypoxia-inducible factor-alpha (HIFα) signaling. Additionally, DE-CMRGs2 significantly inhibited LXR/RXR activation and matrix metalloproteases (Figure 2A and B).

Figure 2 Function enrichment analysis for DE-CMRGs1/DE-CMRGs2. Ingenuity Pathway Analysis (IPA) diagrams of DE-CMRGs1 (A) and DE-CMRGs2 (B).

Identification and Validation of Biomarkers

We took the intersection of DE-CMRGs1 and DE-CMRGs2 to obtain the candidate biomarkers. The Venn diagram revealed that 86 candidate biomarkers were identified as common between DE-CMRGs1 and DE-CMRGs2 (Figure 3A). Among the machine learning models tested, the RF model demonstrated the highest AUC value of 0.953 for distinguishing between samples from patients with CA and NC (Figure 3B). For comparisons between samples from patients with CRC and NC, the GBDT model performed best, achieving an AUC value of 0.935 (Figure S2A). Similarly, the GBDT model also yielded the highest AUC value of 0.929 when comparing samples from patients with CRC and CA (Figure S2B). The top 10 genes from each model were selected and ranked, as shown in Table 2. Overlapping these 30 genes across the models led to the identification of five biomarkers: ZEB1, ABCA1, SLC24A3, CAV1, and FLNA (Figure S2C).

Table 2 Importance Ranks for the Top 10 Model Genes in the Three Optimal Models in CA/NC/CRC

Figure 3 Diagnostic value of the biomarkers in the GSE117606 dataset. (A) Venn diagram showing 86 candidate biomarkers; (B) Diagnostic value of five machine-learning models between CA and NC in GSE117606; (C) Diagnostic value of five biomarkers between CA and NC in GSE117606; (D) Construction of the nomogram based on five biomarkers for clinical use in GSE117606; (E) Statistics for the predictive accuracy of biomarkers on NC, CA, CRC outcomes in both GSE117606 and GSE20916; (F) Expressions of five biomarkers in GSE117606; (G) Expressions of five biomarkers in GSE20916; (H) Diagnostic value of five biomarkers between CA and NC in GSE20916 (ns, not significant, * p <0.05, ** p <0.01, *** p <0.001, **** p <0.0001).

Five diagnostic models were constructed separately using the candidate hub genes and compared based on the AUC values, and the model with the largest AUC value was selected as the best model. The analysis of each individual feature gene’s diagnostic value revealed that ZEB1, ABCA1, SLC24A3, CAV1, and FLNA exhibited strong discriminatory power between CA and NC, as well as CRC and CA, while only CAV1 and FLNA effectively distinguished between CRC and NC (Figure 3C and Figure S2D and E). Additionally, the diagnostic performance of the five-gene signature demonstrated excellent predictive accuracy for NC, CA, and CRC outcomes in both the GSE117606 and external GSE20916 validation sets (Figure 3D). The nomogram based on these five biomarkers, shown in Figure 3E, also exhibited high predictive accuracy for CA, as confirmed by the calibration curve (Figure S2F). The DCA indicated that the nomogram provided clinical benefits across a wide high-risk threshold range (0–1), surpassing the predictive abilities of individual genes such as ZEB1, ABCA1, SLC24A3, CAV1, and FLNA (Figure S2G). Furthermore, the clinical impact curve confirmed the nomogram’s relatively accurate predictive capabilities (Figure S2H). Given the significant predictive value of clinical features, age and gender were incorporated into the five-gene-based nomogram (Figure S3A). This enhanced nomogram demonstrated improved accuracy and clinical benefits compared to individual factors alone (Figure S3B and Figure S3C), with the clinical impact curve further validating these results, consistent with the previous nomogram (Figure S3D).

As illustrated in Figure 3F, all five genes were expressed at lower levels in samples from patients with CA compared to NC and those with CRC. Notably, K-M survival analysis indicated significant differences in survival rates among patients with CRC exhibiting high versus low FLNA expression, where higher FLNA expression correlated with poorer prognosis (Figure S4). In the external validation using the GSE20916 dataset, the expression of four biomarkers (ABCA1, SLC24A3, CAV1, and FLNA) was confirmed, as ZEB1 was not detected. Among these, ABCA1, CAV1, and FLNA exhibited significant differences in expression across samples from patients with CRC, CA, and NC, with FLNA expression aligning with the findings from the GSE117606 dataset (Figure 3G). In the GSE41657 dataset, significant differences were observed in the expressions of CAV1 and SLC24A3 among the samples from patients with CRC, CA, and NC. The expression of CAV1 was consistent with the results of the GSE117606 dataset (Figure S5). Moreover, ROC curve analysis demonstrated that CAV1 and FLNA could effectively distinguish CRC from CA (AUC > 0.7) (Figure 3H, Figure S6A and B). In the GSE41657 dataset, ROC curve analysis demonstrated that CAV1 could be effectively used to distinguish CRC from CA (AUC > 0.7) (Figure S6C–E).

GSEA for Biomarkers

We analyzed the correlation between the biomarkers and the remaining genes to obtain the correlation coefficient. Then, the coefficients of the genes were ranked, and then the ranked results were used to perform GSEA to obtain the biological process of biomarkers enrichment. The GSEA results revealed significant enrichment of various pathways and biological processes associated with the expression levels of CAV1 and FLNA. In the high CAV1 expression group, key pathways such as the Apelin signaling pathway, cyclic guanosine monophosphate-protein kinase G (cGMP-PKG) signaling pathway, adherens junctions, mitogen-activated protein kinases (MAPK) signaling pathway, neuroactive ligand-receptor interactions, proteoglycans in cancer, regulation of the actin cytoskeleton, and vascular smooth muscle contraction were notably enriched. Conversely, in the low CAV1 expression group, neuroactive ligand-receptor interaction were significantly enriched (Figure 4A). For FLNA, the low-expression group showed significant enrichment in ribosomal processes and ribosome biogenesis in eukaryotes. In contrast, the high FLNA expression group exhibited enrichment in pathways such as the Apelin signaling pathway, cGMP-PKG signaling pathway, adherens junctions, MAPK signaling pathway, oxytocin signaling pathway, proteoglycans in cancer, regulation of the actin cytoskeleton, and vascular smooth muscle contraction (Figure 4B). In summary, these biomarkers, particularly CAV1 and FLNA, were predominantly enriched in pathways related to the cGMP-PKG signaling pathway, Apelin signaling pathway, focal adhesion, proteoglycans in cancer, vascular smooth muscle contraction, Salmonella infection, neuroactive ligand-receptor interaction, and regulation of the actin cytoskeleton (Figures 4C–E).

Figure 4 Single gene set enrichment analysis (GSEA) for five biomarkers. (A) CAV1; (B) FLNA; (C) ABCA1; (D) SLC24A3; (E) ZEB1.

Effect of Biomarkers on Immune Heterogeneity

Using ssGSEA to calculate the 24 TME cell score for each sample, the composition of 24 TME cell types in samples from patients with NC, CA, and CRC is illustrated in Figure 5A. Analysis of immune cell proportions revealed significant differences in B cells, macrophages, mast cells, and T cells between samples from patients with CA and NC, as identified by the rank-sum test (Figure 5B). Spearman correlation analysis further highlighted significant positive correlations between the biomarkers and mast cells, as well as eosinophils. Additionally, CAV1, ABCA1, and ZEB1 showed significant positive correlations with both B cells and T cells (cor > 0.2, adj.p < 0.05). In contrast, FLNA exhibited significant negative correlations with macrophages (cor < −0.1, adj. p < 0.05) (Figure 5C).

Figure 5 Immune heterogeneity analysis. (A) Composition of the 24 immune cell types in samples from patients with CA, CRC, and NC; (B) Differences in the proportions of immune cells between different groups, yellow for CA, red for CRC, blue for NC; (C) Correlation network of biomarkers and immune cells, connected lines represent relationships with mutualisms (ns, not significant, * p <0.05, ** p <0.01, *** p <0.001, **** p <0.0001).

Identification of Feature Gene Expression in Different Cell Subsets

Using single-cell RNA-seq to analyzed the expression patterns of biomarkers during cell state transitions and changes in cells during disease states, after data filtering, the features of the scRNA-seq dataset are presented in Figure S7A. A total of 2,000 highly variable genes were selected for cell type identification (Figure S7B). Following normalization, the results of PCA are shown in Figure S7C, from which 30 principal components were selected for further analysis based on a p-value of < 0.05 (Figure S7D). Before clustering, a resolution of 0.06 was determined to be optimal, resulting in well-defined clusters (Figure S7E and Figure S7F). The visualization of the 13 clusters using uniform manifold approximation and projection (UMAP) and t-distributed stochastic neighbor embedding (t-SNE) is depicted in Figure 6A. Heat maps displaying the top five marker genes for each cluster were generated (Figure 6B).

Figure 6 Single-cell RNA-sequencing analysis to identify the different cell subsets in GSE201348. (A) t-Distributed Stochastic Neighbor Embedding (tSNE) plot and Uniform Manifold Approximation and Projection (UMAP) plot for 13 principal components (PCs) identified using principal component analysis (PCA). Each PC captures some biological differences, with the earlier PCs containing more information about these differences; (B) Expression heatmap for the top 5 marker genes among 13 clusters; (C) t-SNE plots labeling cluster 2 as goblet cells based on the expression patterns of FCGBP, MUC2, and CLCA1; (D) t-SNE plots labeling cluster 3 as myeloid cells based on the expression patterns of FTL, TMSB10, and TPT1; (E) t-SNE plots labeling cluster 4 as B cells according to the expression patterns of BLNK, MS4A1, and CD37; (F) A pie chart showing the proportion of 10 cell subsets in samples from patients with NC, CA, and CRC.

Cell types were then identified for each cluster. Cluster 2, characterized by higher expression levels of marker genes such as FCGBP, MUC2, and CLCA1, was identified as the goblet cell cluster (Figure 6C). Cluster 3, with elevated expression of marker genes like FTL, TMSB10, and TPT1, was identified as the myeloid cell cluster (Figure 6D). Cluster 4, showing higher expression of marker genes such as BLNK, MS4A1, and CD37, was identified as the B cell cluster (Figure 6E). Consequently, the 13 clusters were further classified into 10 distinct cell subsets (Table S3). The distribution of cell counts across each sample is shown in Figure 6f, highlighting the presence of epithelial cells, goblet cells, myeloid cells, B cells, fibroblasts, LGR5+ stem cells, lymphocytes, dendritic cells, macrophages, and neurons as the predominant cell types in samples from patients with NC, CA, and CRC.

Dot plots were then employed to visualize the marker genes for each cell type, as depicted in Figure 7A. The distribution of the 10 identified cell subsets is presented in Figure 7B. Additionally, a heat map and a violin plot were used to illustrate the expression levels of the biomarkers within each cell subset, as shown in Figure 7C and D. It is noteworthy that nearly all biomarkers were predominantly expressed in fibroblasts, with FLNA displaying the highest expression levels (Figure 7C and D). Furthermore, the t-SNE cluster map provided a detailed visualization of feature gene expressions across fibroblasts, lymphocytes, goblet cells, B cells, and macrophages (Figure 7E).

Figure 7 Single-cell RNA-sequencing analysis for biomarkers. (A) Dot plots of marker gene expressions for 10 cell subsets, the larger the diameter, the higher the expression; (B) t-SNE plots for 10 cell subsets; (C) Expression heatmap for five biomarkers among 10 cell subsets; (D) Violin plot for feature gene expression in 10 cell subsets; (E) t-SNE cluster map of five biomarkers in different cell subsets, the bluer the color, the higher the expression level.

ceRNA Network and Potential Therapeutic Drugs for Biomarkers

A comprehensive analysis of the five biomarkers identified 57 miRNAs and 21 lncRNAs, which facilitated the construction of a ceRNA network (Figure 8A). The findings indicated that AC138035.1 modulates CAV1 expression via hsa-let-7a-5p or hsa-let-7b-5p, AC240565.2 potentially influences ZEB1 expression through hsa-let-7b-5p, and LINC01001 similarly regulates CAV1 via hsa-let-7b-5p. Additionally, utilizing the DGIdb database, 12 potential drugs were identified targeting the four biomarkers (ABCA1, ZEB1, FLNA, and CAV1). Notably, the analysis revealed a strong binding affinity between ABCA1 and probucol, ZEB1 and salinomycin, as well as FLNA and PTI-125 (Figure 8B).

Figure 8 Exploration of the biomarkers. (A) LncRNA-miRNA-mRNA network relevant to biomarkers; (B) Feature gene-associated targeted drug network using the Drug-Gene Interaction database (DGIdb) database. The pink circle represents the biomarkers, the blue circle represents the targeted drugs, the lines represent the interaction scores, and the color of the lines deepen as the interaction score increases.

Verification of Feature Gene Expression Using qPCR and the HPA Database

The expression of the biomarkers (ZEB1, ABCA1, SLC24A3, CAV1, and FLNA) was validated through qPCR and analysis using the HPA database. qPCR results revealed that the mRNA levels of these genes in samples from patients with CA were significantly lower compared to those with NC (Figure 9), corroborating the findings from the bioinformatics analysis. Additionally, protein levels of the five biomarkers were reduced in samples from patients with CRC relative to those with NC, consistent with the expression patterns observed in the publicly available GSE117606 dataset (Figure 10A–E).

Figure 9 Validation of feature gene expression through the quantitative polymerase chain reaction (qPCR) experiments. * p <0.05, *** p <0.001.

Figure 10 Exploration of the protein expressions of the biomarkers in CRC tissue samples and normal colon samples. (A) ZEB1; (B) ABCA1; (C) SLC24A3; (D) CAV1; (E) FLNA.

Discussion

Early detection of CA is essential for preventing the progression to CRC. Thus, identifying CMRGs linked to both CA and CRC is of substantial clinical importance, as it enables timely interventions. The role of copper metabolism and homeostasis in the development of CA and CRC has garnered increasing interest. Building on the potential relevance of CMRGs in colon cancer,^22,23 a comprehensive bioinformatics analysis was undertaken to identify key genes associated with copper metabolism that contribute to the progression of CA and CRC.

Through this analysis, five biomarkers—ZEB1, ABCA1, SLC24A3, CAV1, and FLNA—were identified. Evaluating the diagnostic value of each gene individually and collectively revealed that ZEB1, ABCA1, SLC24A3, CAV1, and FLNA possess strong discriminatory power between CA and NC.

Zinc finger E-box binding homeobox 1 (ZEB1) plays a pivotal role in the epithelial-mesenchymal transition (EMT) of tumor cells, a process essential for metastasis in colon cancer cell lines.^37,39–41 Silencing ZEB1 has been shown to inhibit the proliferation and metastatic potential of colon cancer cell lines, such as LoVo.⁴² Studies have shown that during the adenoma (CA) stage, ZEB1 expression is generally downregulated, which may be related to maintaining normal epithelial architecture and function, thereby restricting abnormal cell migration and proliferation.⁴³ However, as CA progresses to colorectal cancer (CRC), ZEB1 expression increases significantly.⁴⁴ Elevated ZEB1 expression promotes epithelial–mesenchymal transition (EMT), facilitating tumor cell detachment from the primary site and contributing to cancer metastasis.⁴⁵ ZEB1-mediated regulation of colorectal cancer by cancer-associated fibroblasts (CAFs) has revealed that heterogeneous expression and defects in ZEB1 within these fibroblasts not only compromise the stability of the extracellular matrix (ECM) and the barrier function of CAFs but also collectively contribute to carcinogenesis, suppression of anti-tumor immunity, and promotion of metastasis. This occurs through the upregulation of cytokines mediated by Nuclear Factor kappa-light-chain-enhancer of activated B cells (NF-κB) and the activation of immune checkpoints, underscoring the multifaceted role of ZEB1 in the tumor microenvironment and its significance as a potential therapeutic target.⁴⁶ Copper played a critical role in cellular redox reactions, and its metabolic imbalance could lead to the generation of oxidative stress.⁴⁷ ZEB1 has been shown to repress E-cadherin (CDH1), thereby promoting epithelial–mesenchymal transition (EMT), and may also influence the expression of antioxidant enzymes such as superoxide dismutase (SOD) and glutathione peroxidase (GPX), altering the cellular redox state.⁴⁸ In colorectal cancer (CRC), oxidative stress may contribute to disease pathogenesis by affecting processes such as intestinal epithelial cell proliferation, apoptosis, and inflammation.⁴⁹ Therefore, dysregulated copper metabolism may promote the progression of colorectal adenoma (CA) and CRC through oxidative stress, with ZEB1 potentially serving as a mechanistic link in this process. Our study observed that ZEB1 expression was downregulated in CA samples but upregulated in CRC, reflecting its dual role in cancer progression. However, this phenomenon warrants further investigation. ATP-binding cassette transporter A1 (ABCA1) is a transmembrane protein that mediates the reverse transport of cholesterol from cells to the circulatory system.^50,51 Wu et al proposed that inducing ABCA1 expression could be a novel cancer treatment strategy by reducing intracellular cholesterol levels through reverse cholesterol transport.⁵² In the tumor microenvironment, copper metabolism dysregulation can lead to intracellular changes that may indirectly affect ABCA1-mediated cholesterol transport, thereby influencing tumor cell growth and metastasis.^53,54 Although the direct relationship between ABCA1 and epithelial–mesenchymal transition (EMT) remains unclear, aberrant cholesterol metabolism can alter membrane properties and impact cell–cell adhesion, potentially interacting with EMT processes. Copper-mediated regulation of ABCA1 may serve as a mechanistic link in this context.^55,56 Abnormal ABCA1 expression has been documented in colon cancer tissues and cells, where dysregulated intracellular cholesterol, driven by ABCA1 overexpression in advanced CRC, may contribute to tumor growth and metastasis. Conversely, silencing ABCA1 has been found to promote colon cancer cell proliferation and inhibit apoptosis.^50,57 This study similarly observed downregulation of ABCA1 in patients with CA, which reversed during the progression to CRC, paralleling the pattern seen with ZEB1. This phenomenon supports a potential role for ABCA1 in tumour progression. Caveolin-1 (CAV1), a membrane protein implicated in various cellular processes including endocytosis, extracellular matrix organization, and cholesterol distribution, has been identified as an oncogenic factor in several cancers, including liver, colon, breast, kidney, and lung cancers.^58,59 However, CAV1 expression varies depending on tumor type and stage. While it is downregulated in certain sarcomas and adenocarcinomas, it is significantly upregulated in advanced tumors like squamous cell carcinoma, a finding consistent with our results.⁵⁹ It is important to highlight that there are various copper - transport proteins, such as ATP7A and ATP7B, within cells to maintain copper homeostasis. Given the close relationship between the cell membrane environment housing CAV1 and these transport proteins, it is highly possible that CAV1 interacts with copper - transport proteins.⁶⁰ Research indicates that copper - homeostasis imbalance, either excess or deficiency, can lead to pathological changes like atherosclerosis.⁶¹ These links imply a potentially complex connection between CAV1, copper - transport proteins, and copper homeostasis in tumor - related physiological and pathological processes. In colorectal adenoma (CA), CAV1 expression is usually low. Low CAV1 expression may limit the activation of cell - proliferation and migration - promoting signaling pathways, keeping the cell relatively stable.⁶² Our study shows that low CAV1 expression is significantly related to neural - activity ligand - receptor interaction and the ribosome pathway. Neural - activity ligands like neurotransmitters may have their receptors’ expression and function on the cell surface influenced by copper - ion levels.⁶³ When CA progresses to colorectal cancer (CRC), CAV1 expression increases. High CAV1 expression may boost endocytosis, affecting the number and distribution of cell - surface receptors and thus cell response to growth - factor - like signaling molecules.⁶⁴ Notably, Friedrich et al reported that the absence of CAV1, especially when combined with certain tumor suppressors and adenomatous polyposis in mouse models, may increase the risk of CRC.^59,65 In the study by Li et al, HSPA8 was found to be transcriptionally upregulated in BRAF V600E mutant CRC, promoting the degradation of CAV1 and the release of β-catenin into the nucleus, thereby activating the Wnt/β-catenin pathway and driving metastasis and progression in BRAF V600E CRC.⁶⁶ However, the expression of CAV1 appears to be context-dependent. In our study, we observed that CAV1 expression was downregulated in colon adenomas (CA) but upregulated in CRC. Filamin A (FLNA), the most abundant and broadly distributed member of the actin-binding protein family, has diverse functions.^67,68 As essential cofactors for numerous enzymes, copper ions participate in multiple intracellular metabolic pathways and may indirectly modulate the expression and function of Filamin A (FLNA).^54,69 Our study reveals a significant association between ribosomal biogenesis pathways and FLNA underexpression. Copper-dependent enzymes, such as copper-zinc superoxide dismutase (SOD1), may contribute to rRNA stability or processing, thereby influencing ribosomal assembly.⁷⁰ In colorectal cancer (CRC), dysregulated copper metabolism could potentially drive tumor progression through ribosomal pathway activation.⁷¹ Furthermore, within the tumor microenvironment, altered copper ion concentrations may indirectly affect FLNA’s role in epithelial-mesenchymal transition (EMT) by modulating cytoskeleton-associated signaling pathways.⁷² In the context of CRC, the clinical expression and biological role of FLNA remain controversial. Tian et al found that FLNA is downregulated in patients with CA and is closely associated with CRC development and progression, suggesting it as a prognostic marker.⁷³ Conversely, other research indicates that FLNA overexpression and phosphorylation at ser2152 are linked to snail-induced EMT and cell adhesion, proposing FLNA as a potential target for combined therapies in CRC progression.⁷⁴ This could explain the observed upregulation of FLNA in CA and its downregulation in CRC samples in this study. Cheng demonstrated that FLNA directly regulates metastasis and epithelial-mesenchymal transition (EMT) in chemotherapy-resistant CRC cells, suggesting its potential as a reliable CRC marker. This is consistent with our findings, which identified FLNA as a primary marker for CRC diagnosis.⁷⁵ Additionally, SLC24A3, a sodium-potassium-calcium exchanger, is associated with poor prognosis when expressed at high levels, indicating its utility as a prognostic marker, particularly in cervical cancer.^76,77 This study is the first to suggest a potential role for SLC24A3 in the initiation and progression of colon cancer. However, further research is required to elucidate the roles and underlying molecular mechanisms of these five characteristic genes in CA and CRC.

To further elucidate the relationship between characteristic genes and key immune cells in CA, immune cell infiltration was analyzed. Previous studies, such as those by Mo et al, have shown that the downregulation of the ITGA4 gene is positively correlated with increased proportions of mast cells and eosinophils in CA.⁷⁸ The protective and antitumor roles of mast cells in early-stage intestinal tumor models have also been demonstrated.⁷⁹ Consistent with these results, the biomarkers identified in this study exhibited significant positive correlations with the infiltration levels of these immune cells. Additionally, the association between a higher proportion of M0 macrophages, low FLNA expression, and poor prognosis in CRC, as documented in previous literature,⁸⁰ was confirmed by the current study in CA samples. These results underscore the close relationship between characteristic genes and immune cell infiltration in both CA and CRC. Subsequent scRNA-seq analysis revealed that these biomarkers were expressed across various cell types, including fibroblasts, lymphocytes, goblet cells, B cells, and macrophages. Notably, fibroblasts exhibited the highest expression of nearly all characteristic genes, with FLNA showing particularly elevated levels. The critical role of fibroblasts in tumor initiation and progression, where fibroblast-related characteristics serve as prognostic markers, has been well-documented.⁸¹ Hiroki Kobayashi et al further highlighted the importance of cancer-associated fibroblasts in CRC progression and their link to poor patient prognosis.⁸² These findings reinforce the diagnostic value of the selected characteristic genes in CRC.

To delve deeper into the pathways and functions associated with the identified characteristic genes, GSEA was conducted. The analysis revealed significant enrichment of these genes in several pathways, including the cGMP-PKG signaling pathway, Apelin signaling pathway, focal adhesion, and proteoglycans in cancer. Therapeutic activation of the cGMP-PKG pathway has been proposed as a promising strategy for colon cancer prevention and treatment, with PKG playing a central role in intestinal homeostasis and tumor suppression by inhibiting the β-catenin, T-cell factor (TCF), and SOX9 signaling pathways, thereby blocking cell proliferation and tumor angiogenesis.⁸³ Additionally, increasing evidence supports the anti-tumor signaling capabilities of PKG1 in colon cancer cell lines.⁸⁴ The Apelin pathway, which activates G-protein coupled receptors and Apelin receptors, also exerts various biological functions.^85,86 Notably, Picault et al found that the Apelin gene was upregulated in approximately half of the CA cases, suggesting that the Apelin signaling pathway may contribute to colon cancer progression by activating anti-apoptotic pathways.⁸⁷ The focal adhesion signaling network, composed of integrins, growth factor receptors, and other survival-promoting molecules, represents another potential target for cancer treatment.⁸⁸ Moreover, studies have shown that syndecan-2, a cell surface proteoglycan, is upregulated in CRC and promotes cell migration, while the expression of the antineoplastic molecules syndecan-1 and syndecan-4 is reduced.⁸⁹ These results suggest that the dysregulation of these biological processes is closely linked to CRC onset and development, highlighting the potential of the five biomarkers as indicators of pathways associated with CRC progression in patients with CA. The incorporation of these genes into diagnostic strategies and treatment plans holds significantly to enhance our ability to manage CRC. Their expression levels in specific cell types, notably fibroblasts and immune cells, suggest that they may also significantly influence the tumour microenvironment and impact response to treatment.

Moreover, an upstream regulatory network involving lncRNA-miRNA-mRNA interactions was constructed for the characteristic genes linked to disease progression. Vahid et al proposed that MALAT1, regulated by the lncRNA-p53 network, could serve as an early detection biomarker for CRC.⁹⁰ In this study, MALAT1 was found to competitively bind to miR-141-3p, thereby influencing the expression of SLC24A3. Additionally, Mezgebe et al demonstrated that miRNA let-7c, which interacts with CAV1 and ZEB1 in our study, was significantly upregulated in adenomas compared to healthy samples.⁹¹ This finding may explain the downregulation of CAV1 and ZEB1 in CA samples. The constructed network offers a foundation for further exploration of the underlying mechanisms governing these genes in CA and CRC.

Among the drugs predicted to target the characteristic genes involved in disease progression, probucol, which targets ABCA1, has been shown to effectively inhibit ABCA1-mediated cholesterol efflux,⁹² a process implicated in the development of various tumors. CAV1 demonstrates a strong binding affinity for testosterone, and studies have shown that both endogenous and exogenous testosterone can exacerbate azoxymethane (AOM) or dextran sodium sulfate (DSS)-induced colitis and its carcinogenic effects.⁹³ Another potential drug, salinomycin, targets ZEB1 and has shown promise in eradicating cancer stem cells (CSC) and treating drug-resistant cancer cells, suggesting its potential as an effective chemotherapy agent.⁹⁴ Wang et al found that salinomycin significantly reduced tumor growth and CSC-related Wnt target gene expression, including LGR5, in colorectal tumor xenografts and APCmin/+ transgenic mice.⁹⁵ Therefore, the predicted drugs highlighted above hold considerable promise as potential therapeutic candidates for the treatment and prevention of CA and CRC.

While our study provides valuable insights, it is important to acknowledge its limitations. The current sample may not fully capture the heterogeneity of all colon adenoma (CA) patients. While this study provides valuable insights into the roles of copper metabolism and the immune microenvironment in CA, several limitations should be acknowledged. The current sample cohort may not fully capture the heterogeneity of CA patients. The qPCR validation was performed with a limited sample size and lacks long-term clinical follow-up. Subsequent studies should expand the sample size to include individuals of diverse ethnicities, ages, genders, and clinical stages, coupled with longitudinal monitoring of gene expression patterns and disease progression to enhance the generalizability and reliability of the findings. Furthermore, the stability of signature gene expression may vary across different disease stages, necessitating longitudinal studies to validate their clinical utility. We plan to systematically analyze existing datasets, incorporating clinical staging, family history, and other variables to construct a nomogram for evaluating associations and interactions with the five key genes, thereby improving predictive accuracy. Addressing these limitations is essential for translating these research findings into clinical applications.

To enhance the scientific and translational value of the study, animal models and cell culture experiments should be conducted in the future to validate the biological functions of these signature genes by knocking out or overexpressing them, thereby confirming their roles in the pathogenesis of CA. Secondly, the dynamic trends of cancer characteristics were effectively captured by monitoring the biomarkers through longitudinal studies using the Smoothing Cubic B-Spline Model,⁹⁶ which combines multi-omics data and reinforcement learning state behaviour reward state Behaviour (SARSA) hybrid cancer prediction model, which lays a solid foundation for subsequent bioinformatics analyses.^97–99 Combined with multi-centre and large-sample clinical studies, we will assess its clinical application value in different stages of the disease and enhance the generality of the findings.

Through the efforts in the above directions, we expect to overcome the limitations of current studies, further reveal the molecular mechanisms of CA, promote clinical translation, and open up new avenues for precision diagnosis and treatment of CA.

Conclusions

In conclusion, this study successfully identified and validated Cu metabolism-related biomarkers associated with CA and CRC through comprehensive bioinformatics approaches. The research also examined the correlation between the biomarkers and immune infiltration, alongside predicting potential targeted drugs for CA treatment. Additionally, a ceRNA regulatory network was constructed to elucidate lncRNA-miRNA-mRNA interactions, offering valuable insights for future mechanistic studies. Overall, this study lays a solid foundation for understanding the molecular mechanisms underlying these diseases.