Matrix data from GSE100480 (19), GSE49355 (20), GSE81558 (21), and GSE41258 (22) were downloaded from the Gene Expression Omnibus (GEO) database. Expression data of GSE100480 and GSE81558 were normalized using robust multichip average (RMA) method; thus, log₂ transformation not required. GSE49355 was normalized using MAS5 method and GSE41258 using PLIER method before being uploaded; log₂ transformation were performed to scale these data. Phenodata were investigated, and expression data of liver metastasis (LM) and primary colorectal cancer (PC) were extracted for later analysis. Quantile normalization were applied using R package preprocessCore to ensure the data have the same distribution (23).

Least squares, empirical Bayes, and t-test methods based on R package limma were used to analyze the DEGs between LM and PC in four GSEs separately (24). Probes representing multiple genes and duplicated genes were omitted. A p < 0.05 and | log₂fold change (FC) | >1 were defined as the threshold for DEGs screening. Robust rank aggregation method from R package RobustRankAggreg was executed to integrate DEGs from four GSEs (25). Adjusted p < 0.05 and | log₂FC | > 1 were set as the criteria to filter statistically significant DEGs.

Gene set enrichment analysis (GSEA) was performed in R package clusterProfiler (26, 27), with genes metric ranked according to average log₂FC and Kyoto Encyclopedia of Genes and Genomes (KEGG) gene sets from the Molecular Signatures Database (MSigDB) as input. p < 0.05 and adjusted p < 0.25 were considered as significant.

KEGG pathway and Gene Ontology (GO) of integrated up- and downregulated genes were enriched and annotated by the Database for Annotation, Visualization and Integrated Discovery (DAVID), respectively, including biological process (BP), cellular component (CC), and molecular function (MF) (28, 29). Top terms with p < 0.05 were deemed as significant and visualized with R package ggplot2 (30).

Relationships among integrated DEGs were evaluated by STRING (31); interactions with combined score >0.4 were exported to Cytoscape (version 3.8.2, https://cytoscape.org). Nodes’ scores were calculated by plug-in cytoHubba using 12 methods, and the top 50 genes of each method were kept. Genes existing at least 10 out of 12 methods were selected as candidate hub genes.

The Cancer Genome Atlas (TCGA) colorectal cancer cohort were divided into two groups according to median gene expression level. Gene Expression Profiling Interactive Analysis 2.0 (GEPIA2) was used to evaluate the difference in disease-free survival (DFS) between the groups (32). Mantel–Cox test value <0.05 was set to determine hub genes associated with poor prognosis of colorectal cancer. Comparisons of hub genes’ expression between LM and PC from GSE41258 were analyzed using Wilcoxon rank sum test with continuity correction; p < 0.05 were considered as statistically significant.

The online website oPOSSUM 3.0 was used to predict the transcription factor (TF) of the hub genes (33). Overrepresented conserved TF binding sites were detected based on the criteria that the conservation cutoff was 0.4. The amount of upstream/downstream sequence was 5,000/5,000; the species was Homo sapiens, with Z-score value >10. Pearson’s correlation coefficient were used to the determine the relationship between transcription factors and genes; p < 0.05 was considered as statistically significant.

Protein expression level in different stages of CRC was curated from the Clinical Proteomic Tumor Analysis Consortium (CPTAC) by UALCAN. The immunohistochemical images of the hub genes in advance and early stages were also identified using Human Protein Atlas (34, 35). The Drug Gene Interaction Database (DGIdb) was used to predict the potential drugs targeting hub genes (36). Drugs with combined value of query score and interaction score >10 were selected for docking. Homologous structures of gene-encoded protein were downloaded from the Protein Data Bank (PDB) and three-dimension structures of drug from PubChem (37, 38). Autodock (version 4.2.6) software was used to preprocess and define the proteins and drugs to receptors and ligands, respectively (39). The docking grid box was set to envelop the whole receptors; docking parameters was set as genetic algorithm with short maximum number of evaluations. The docking results were ranked by energy, and the first model was exported to Pymol for visualization (40).

HCT116 and LoVo cells were harvested and resuspended with Dulbecco’s modified Eagle’s medium (DMEM) at a concentration of 3 × 10⁶ cells/ml. Then, 200 μl cell suspension with 200 μM ADH-1 (Cat. No. HY-13541, MedChemExpress, USA) or epigallocatechin (Cat. No. HY-N0225, MedChemExpress, USA) was added into the top chambers, and the bottom chambers were filled with DMEM with 10% fetal bovine serum (FBS). Cells were allowed to migrate for 24 h. Non-migrated cells were erased with a cotton swab, and migrated cells were fixed with 4% paraformaldehyde for 20 min and stained with crystal violet.

LO2, Huh7, and SK-Hep1 cells were digested with trypsin to obtain a single cell suspension. Approximately 500 cells were seeded onto a 48-well plate with different concentrations (0, 50, 100, 200, and 400 μM) of ADH-1 or epigallocatechin and incubated for 10 days. When the colony was visible to the naked eye, the colonies were carefully washed twice with phosphate-buffered saline (PBS). Then, the colonies were fixed with 4% paraformaldehyde for 20 min and stained using crystal violet.

LO2, Huh7, and SK-Hep1 cells were digested with trypsin to obtain a single cell suspension. Approximately 2,000 cells were seeded onto a 96-well plate with different concentrations (0,12.5, 25, 50, 100, 200, 400, and 800 μM) of ADH-1 or epigallocatechin and incubated for 72 h. Ten microliters of CCK8 reagent was added and treated for 3 h; then, the absorbance was measured at 450 nm.

R 3.6.3 (https://www.R-project.org/) and GraphPad Prism 5 were used for statistical analysis. In vitro data were expressed as mean ± SEM, Students’ t-test was used for calculation of p values. *** indicates p < 0.0001.

This study was conducted based on the process described in Figure 1 . Briefly, we downloaded expression matrix of four GSE series containing liver metastasis (LM) and primary colorectal cancer (PC) from GEO and employed preprocess and quality control to obtain analysis-ready data. GSEA was performed to determine concordant differences between LM and PC. Next, we applied a robust method to integrate and identify common DEGs between LM and PC. To further investigate the characteristic of LM, GO and KEGG annotation were performed. Meanwhile, we constructed a PPI network to analyze the relationships among DEGs. Along with survival analysis from TCGA colorectal patients, we identified four hub genes upexpressed in LM compared with PC, which were associated with poor prognosis of patients. Furthermore, we predicted transcription factor (TF) of hub genes and estimated their correlation. Simultaneously, potential drugs targeting hub genes were predicted, and protein–drug interaction was evaluated. Lastly, in vitro experiments were conducted to analyze the antitumor abilities of the two predicted drugs.

To diminish the system bias in microarray data and make sure that the difference was biologically significant, quantile normalization method was used in the selected sample of four GSEs ( Supplementary Figures S1A–D ). There were 8 LM and 13 PC in GSE100480, 19 LM and 20 PC in GSE49355, 19 LM and 23 PC in GSE81558, and 47 LM and 186 PC in GSE41258; a total of 93 liver metastasis tumor and 242 primary tumor expression data were included in this study. Limma identified 336 significantly upexpressed and 83 downexpressed genes in GSE100480 ( Supplementary Figure S2A ), 171 upexpressed and 127 downexpressed genes in GSE49355 ( Supplementary Figure S2B ), 189 upexpressed and 69 downexpressed genes in GSE81558 ( Supplementary Figure S2C ), and 88 upexpressed and 102 downexpressed genes in GSE41258 ( Supplementary Figure S2D ). Heatmaps demonstrating expression of DEGs were also shown ( Supplementary Figures S3A–D ). Robust rank aggregation method was executed to integrate DEGs from the four GSE series; eventually, 138 overlapping genes including 108 significantly upexpressed and 30 downexpressed genes were obtained ( Figure 2 and Supplementary Table 1 ).

We performed GSEA to investigate whether a priori defined gene sets relevant to cancer metastasis were significantly different between LM and PC. A total of 22,750 genes were involved. Results suggested that genes comprising the Wnt signaling pathway ( Figure 3A ) and TGF-β ( Figure 3B ), two canonical pathways widely accepted for regulating metastasis of cancer, were not statistically significant. However, we observed significant enrichment in cell adhesion molecules ( Figure 3C ) and PPAR signaling pathway ( Figure 3D ), which may participate in the regulatory role of liver metastasis.

To explain the biological difference between LM and PC groups, functional enrichment of integrated DEGs was also performed. According to the results, upexpressed genes were significantly enriched in the negative regulation of endopeptidase activity of BP, blood microparticle of CC, and serine-type endopeptidase inhibitor activity of MF ( Supplementary Figure S4A ). Besides, enriched KEGG pathway of upexpressed genes were significant in complement and coagulation cascades ( Supplementary Figure S4B ). In addition, the downexpressed genes were mainly associated with the collagen catabolic process of BP, proteinaceous extracellular matrix of CC, and metalloendopeptidase activity of MF ( Supplementary Figure S4C ). Meanwhile, KEGG enrichment of downexpressed genes were significantly participated in pancreatic secretion ( Supplementary Figure S4D ).

The protein–protein interaction network of DEGs in LM was constructed via STRING. Except for combined score <0.4, a total of 138 nodes and 1,388 edges were established in Cytoscape ( Figure 4A ). A node represents a gene encoded protein, and an edge indicates mutual interaction. Plug-in cytoHubba was used to calculate nodes’ scores with 12 methods, and 35 nodes with 506 edges were selected as candidate hub genes ( Figure 4B and Supplementary Table S2 ), which were all remarkably upexpressed in the LM group ( Figures 5A–D ). We next looked into DFS information of TCGA CRC cohort using GEPIA2. High expression of ALB ( Figure 5A ), APOE ( Figure 5B ), CDH2 ( Figure 5C ), and ORM1 ( Figure 5D ) were identified as biomarkers significantly related to poor prognosis of CRC patients. GO and KEGG analysis of these four genes indicated that they were responsible for extracellular exosome and protein binding ( Supplementary Figure S5 ). The above three characteristics illustrated that these four hub genes contribute to LM and severely threat patients’ health.

To further investigate the underlying mechanism that regulates hub genes, potential transcription factors (TFs) were predicted by oPOSSUM 3.0 ( Supplementary Table S3 ). Based on prescribed criteria, we obtained three TFs, namely, ESR2, FOXO3, and SRY. The estimated binding site of TFs and genes was established, respectively ( Figure 6A and Supplementary Table S4 ). We then examined whether an association existed between TFs and hub genes in samples recruited in this study ( Figure 6B ). There was indeed a positive correlation between ALB and FOXO3 expression and that of CDH2 and FOXO3 expression ( Figure 6C ). We also found a negative correlation between APOE and ESR2 expression and ORM1 and SRY expression ( Figure 6D ). Interestingly, the four hub genes were strongly positively correlated with each other (R > 0.5), implying that they might have synergistic effect to promote LMCRC.

Four genes, namely, ALB, APOE, CDH2, and ORM1, were considered as hub genes and highly expressed in LM. We examined protein expression of these genes and found that higher expression accompanied with advanced stage ( Figures 7A–D ). It was consistent with our previous findings. Thus, hub genes might serve as potential therapeutic targets. In this context, we used DGIdb online database to search drug–gene interactions. A total of 58 drugs targeting hub genes were obtained ( Figure 8A and Supplementary Table S5 ). Among them, four drugs, namely, ADH-1, epigallocatechin, CHEMBL1945287, and cochinchinenin C, receiving higher scores than preset criteria were selected. The molecular structures of drugs and proteins downloaded from PubChem (CID: 9916058, 72277, 18877472, and 23634528) and PDB (6RJV and 2QVI) were demonstrated ( Supplementary Figures S6A–F ) separately. Docking results showed that interactions existed between ALB and epigallocatechin ( Figure 8B ) with binding energy of −3.89, CDH2 and ADH-1 ( Figure 8C ) with binding energy of −5.04, ALB and cochinchinenin C ( Supplementary Figure S7A ) with binding energy of −2.32, and ALB and CHEMBL194528 ( Supplementary Figure S7B ) with binding energy of −3.9, respectively. These drugs were potentially prospective agents to prevent and treat LM by interfering with hub genes ALB and CDH2.

To further confirm our findings, in vitro experiments were conducted with two available drugs. First, the capacity of ADH-1 and epigallocatechin on colorectal cancer metastasis was examined. The transwell assay results showed that both ADH-1 and epigallocatechin significantly inhibited the migration of HCT116 and LoVo cells ( Figures 9A, B ). Meanwhile, we also evaluated whether ADH-1 and epigallocatechin had an effect on the proliferation of normal liver cell LO2 and tumor cells such as Huh7 and SK-Hep1. As we expected, colony formation and CCK8 indicated that ADH-1 did not affect the proliferation of normal liver cells and tumor cells ( Figures 9C, D ), while epigallocatechin reduced liver cells’ proliferation abilities ( Figures 9E, F ). These results supported that both potential drugs can prevent colorectal cancer cells metastasis.