Identification of colorectal cancer–associated target proteins

Colorectal cancer–associated target genes were systematically collected from two publicly available disease-related databases, namely DisGeNET ( https://www.disgenet.org/) and GeneCards ( https://www.genecards.org/), accessed on December 21, 2025. ^{11, 12} The search strategy employed the keywords “colorectal cancer” and “colorectal carcinoma” in both databases. All identified gene entries were compiled, after which duplicate records were removed. The curated gene lists from both sources were subsequently integrated to generate a consolidated dataset representing colorectal cancer–related target proteins.

Prediction of vitamin D3-related target proteins

The SMILES string of VD3 was retrieved from the PubChem Database, then submitted independently to SwissTargetPrediction ( https://www.swisstargetprediction.ch/) and the Similarity Ensemble Approach (SEA) platform ( https://sea.bkslab.org/), with both analyses restricted to Homo sapiens to ensure species-specific relevance. Protein targets predicted by the two platforms were combined into a single dataset, followed by the removal of duplicate entries.

Construction of the protein–protein interaction network

Potential therapeutic targets of VD3 in CRC were identified by intersecting peptide-associated proteins with breast cancer–related targets. The overlap between the two datasets was determined using the InteractiVenn web tool ( http://www.interactivenn.net/). The intersecting genes were considered candidate targets involved in the pharmacological action of the peptide. These overlapping targets were subsequently imported into the STRING database (version 12.0; https://string-db.org/) to construct a protein–protein interaction (PPI) network. The analysis employed the “Multiple Proteins” option, restricted the organism to Homo sapiens, and applied a highest confidence score of 0.900. The generated interaction network was exported in TSV format and visualized using Cytoscape software (version 3.10.0). Within the PPI network, node size and color intensity reflected degree centrality, while edge thickness represented the combined interaction confidence score. ¹³

Network topological properties and hub gene identification

Topological characteristics of the PPI network were quantitatively evaluated using the CytoHubba plugin implemented in Cytoscape. Hub genes were identified based on the degree centrality algorithm, and the ten nodes with the highest Maximal Clique Centrality (MCC) were selected as key hub targets. ¹⁴ Overall centrality (OC) and clustering analyses were conducted in RStudio using the FactoMineR, factoextra, ggrepel, and ggplot2 packages.

Gene ontology and KEGG pathway enrichment analysis

To explore the biological processes and signaling pathways associated with the pharmacological activity of VD3 in CRC, enrichment analyses were performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID 2022; https://david.ncifcrf.gov/). ¹⁵ The overlapping target genes derived from the PPI analysis were uploaded to the platform, with the species parameter set to Homo sapiens.

Gene Ontology (GO) enrichment analysis was conducted on the Molecular Function (MF), Cellular Component (CC), and Biological Process (BP) categories, while Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was performed to identify relevant signaling pathways. A significance threshold of p < 0.05 was applied to both analyses. Enriched GO terms and KEGG pathways were ranked based on ascending p-values, and the top ten entries were selected for visualization. Graphical representations, including bubble plots and bar charts, were generated using the SRplot online bioinformatics platform ( https://www.bioinformatics.com.cn/), illustrating enrichment significance, gene counts, and pathway involvement related to VD3 in CRC pathophysiology.

Molecular docking analysis

Molecular docking was performed using MOE 2022.02 to evaluate the binding affinity of VD3 toward selected CRC–associated target proteins. Protein crystal structures were retrieved from the Protein Data Bank, followed by removal of co-crystallized ligands, water molecules, and heteroatoms, and energy minimization using the Amber10:EHT force field. Binding sites were predicted using the Site Finder Feature in MOE 2022.02. Docking was conducted using the triangle matcher algorithm with the London dG scoring function. Docked complexes were ranked based on binding free energy (kcal/mol) and refined through force-field minimization. Two-dimensional and three-dimensional visualizations were generated in MOE to characterize key residue interactions at the ligand–protein interface. ¹⁶

Molecular dynamics simulation

Molecular dynamics simulations were performed employing YASARA Dynamics version 4.3.13. Protein–peptide complexes selected for analysis were loaded into the platform through the “Set Target” function and the “Macro & Movie” module available in the Options menu. Simulations were carried out under conditions approximating physiological homeostasis, including a temperature of 310 K, pressure of 1 atm, pH 7.4, and an ionic strength corresponding to 0.9% NaCl. The AMBER14 force field was utilized, and each system was solvated within a cubic periodic boundary box with a 10 Å margin from the protein surface. The duration of each simulation was defined as 50,000 ps (50 ns). Following completion of the trajectories, structural stability was assessed by calculating the root mean square deviation (RMSD) of Cα atoms using the MD_analysis macro. Graphical visualization and data plotting of the RMSD profiles were subsequently conducted using ORIGINPro 2024 (OriginLab, Massachusetts, United States). ¹⁶

A total of 20,115 CRC–related target proteins were identified, whereas 111 target proteins of VD3 were identified from databases. Cross-analysis of the two datasets revealed 110 overlapping targets, representing the shared molecular targets between VD3 and CRC ( Figure 1A).

A total of 110 overlapping target proteins shared between VD3 and CRC were used to construct the PPI network. Interaction data were obtained from the STRING database ( Homo sapiens, confidence score ≥ 0.900) and visualized using Cytoscape software (version 3.10.0). After removal of isolated nodes lacking defined interactions, the final PPI network comprised a densely connected interaction map, reflecting extensive functional associations among the VD3–CRC targets.

As illustrated in Figure 1B, nodes represent target proteins and edges denote predicted or experimentally supported protein–protein interactions. Topological analysis of the PPI network was performed using the cytoHubba plugin in Cytoscape. The MCC algorithm was applied to identify highly influential nodes within the network, and the top 10 hub genes were selected based on their MCC scores ( Figure 2). The identified hub targets included CDK4, AURKA, TOP2A, CDC45, MELK, SRD5A1, CYP17A1, CYP19A1, ESR1, and AR, all of which exhibited high connectivity and central positioning within the VD3–CRC interaction network. Detailed topological parameters of these hub genes, including MCC, Degree, Betweenness, and Closeness centralities, are summarized in Table 1. Among these targets, CDK4 and CYP19A1 displayed the highest degree and MCC values, indicating dominant regulatory roles within the network.

In the BP category, the most significantly enriched terms included signal transduction, G protein–coupled receptor signaling pathway, lipid metabolic process, and positive regulation of cell population proliferation. Additional enriched BP terms comprised steroid metabolic process, regulation of gene expression, response to xenobiotic stimulus, inflammatory response, and positive regulation of cytosolic calcium ion concentration. These processes exhibited high gene counts and low p-values, indicating robust enrichment among the VD3–CRC overlapping targets ( Figure 3A).

Analysis of the CC category revealed that target genes were predominantly associated with the membrane and plasma membrane, followed by the endoplasmic reticulum, receptor complex, and synapse-related structures, including postsynaptic membrane and neuron projection. Other enriched CC terms included the perinuclear region of cytoplasm, dendrite, and cell body fiber, reflecting diverse subcellular localizations of the identified targets ( Figure 3A).

Within the MF category, protein binding emerged as the most significantly enriched term, followed by metal ion binding, nucleotide binding, and ATP binding. Enrichment was also observed for G protein–coupled receptor activity, kinase activity, transferase activity, and specific subclasses such as protein serine/threonine kinase activity and protein serine kinase activity. These MF terms demonstrated notable gene ratios and statistical significance, indicating functional convergence among the VD3–CRC target proteins ( Figure 3A).

The integrated gene–term network visualization further illustrated the relationships between individual target genes and enriched GO categories, highlighting the distribution of targets across multiple biological processes, molecular functions, and cellular compartments. The size of nodes reflected gene counts, while color gradients corresponded to the level of statistical significance, providing an overview of functional enrichment patterns within the VD3–CRC target network ( Figure 3B).

KEGG pathway mapping was performed to visualize the distribution of VD3–related CRC target genes within canonical CRC signaling pathways ( Figure 4). The mapped targets were primarily localized to pathway branches related to cell proliferation and survival, with prominent involvement of AREG, EGFR, and MEK within the ErbB–MAPK signaling axis. These components were positioned at receptor-level and downstream kinase nodes, indicating coverage of multiple signaling tiers within CRC-related pathways.

In addition, VD3–CRC targets were mapped to interconnected signaling modules, including the MAPK cascade and downstream transcriptional regulators associated with proliferative responses. Overall, KEGG mapping highlights that VD3–CRC targets are integrated into key proliferation-related branches of colorectal carcinogenesis rather than being confined to isolated molecular nodes.

Molecular docking was performed to assess the interactions between VD3 and ten hub target proteins associated with CRC, including CDK4, AURKA, CYP17A1, TOP2A, CDC45, SRD5A1, MELK, CYP19A1, AR, and ESR1 ( Table 2). All VD3–protein complexes exhibited negative binding free energy values, ranging from −5.72 to −8.74 kcal/mol. The strongest binding affinities were observed for CYP19A1 (−8.74 kcal/mol), SRD5A1 (−8.23 kcal/mol), and CDC45 (−7.98 kcal/mol). RMSD values for all docked complexes were below 2.0 Å (0.96–1.94 Å), indicating stable and reliable docking conformations. Two-dimensional and three-dimensional visualizations were generated to illustrate the binding modes of VD3 within the target protein pockets ( Figure 5).

Molecular dynamics simulations over 50 ns were performed to assess the stability of vitamin D complexes with CYP19A1, AR, AURKA, CDC45, CDK4, CYP17A1, ESR1, SRD5A1, TOP2A, and MELK ( Figure 6). Backbone RMSD analysis showed initial equilibration followed by stable trajectories for all complexes. AR and CDC45 exhibited the lowest RMSD values (~1.2–1.5 Å), indicating high structural stability, whereas TOP2A, MELK, and ESR1 displayed higher fluctuations, with RMSD values ranging between ~2 and 4.8 Å.