Two comprehensive transcriptome datasets were employed: 1) the microarray dataset (GSE8671) (Sabates-Bellver et al., 2007) which compiled on Affymetrix Human Genome U133 Plus 2.0 Array platform (Affymetrix Inc.,Santa Clara, CA, United States) with 32 colorectal adenoma samples and 32 matched normal tissue samples obtained from NCBI-Gene Expression Omnibus (GEO) (Barrett et al., 2013), and 2) the RNA-seq COAD and READ datasets (normalized as FPKM) from The Cancer Genome Atlas (TCGA) consisting of 644 tumor tissue samples, and 51 normal tissue samples (Tomczak et al., 2015). TCGA-COAD consisted of 478 primary tissue samples and 41 normal colon tissue samples. TCGA-READ consisted of 166 primary tissue samples and 10 normal rectal tissue samples.

BioGRID database (v.3.5.167) (Chatr-aryamontri et al., 2017) was used, which contains 35,688 physically and experimentally detected PPIs among 8,570 human proteins for the human protein interactome. Filtering the interactome dataset with proteins encoded by the genes represented in the transcriptome datasets resulted in a network of 34,603 PPIs among 8,322 proteins considering the TCGA dataset and 32,259 PPIs among 7,951 proteins considering the GSE8671 dataset.

The differential interactome algorithm was recruited to identify the differential PPIs (dPPIs) between the tumor (or adenoma) phenotype and the normal phenotype, considering the relative observation frequencies (q-value) of each PPI as previously described (Ayyildiz et al., 2017; Gulfidan et al., 2020). Briefly, genes were grouped into three different levels as −1, 0, 1 according to their expression levels within each sample. When the expression level of a gene was lower than the average expression value, it was labeled as “−1’; when the expression level of a gene was higher than the average, it was labeled as “1”; otherwise, it was labeled as “0”. The false discovery rate was 0.05. According to the three-level expression categorization, there were nine possible gene expression states (i.e., [0 0], [0 1], [0 −1], [1 0], [1 1], [1 −1], [−1 0], [−1 1], [−1 −1]) for each interacting protein pair. The number of times the conditions occurred in the normal group (N₀) and the number of times they occurred in the tumor group (N₁) were calculated. Taking into account the imbalance between the sample sizes of the groups, the count parameters were normalized considering the total sizes of the normal (N_N) and tumor (N_T) groups (the maximum possible numbers of N₀ and N₁, respectively), and the normalized observation frequencies in each group were obtained. The q value (the estimate of the probability of that state occurring in the tumor state) was calculated as follows: q =   N 1 N T N 0 N N + N 1 N T

PPIs with a normalized observation frequency in either the normal or tumor phenotype greater than 20% and q-values of less than 0.10 (significantly repressed in the tumor phenotype) or greater than 0.90 (significantly activated in the tumor phenotype) were considered dPPIs. Differential interacting proteins (DIPs) were defined as proteins that exhibit significant alterations in their interaction patterns during the transition from normal phenotype to tumor state (Gulfidan et al., 2020), and categorized into two groups depending on their interaction patterns: 1) DIPs with repressed interactions under tumor conditions, and 2) DIPs with activated interactions under tumor conditions. Modules were constructed around DIPs using their dPPIs and visualized using Cytoscape 3.7.2 (Shannon et al., 2003).

Gene set over-representation analyses were carried out using the ConsensusPathDB (ver.34) (Kamburov et al., 2013). Preferred data sources for metabolic and signalling pathways were KEGG (ver.88.0) (Kanehisa et al., 2017), Reactome (ver.64) (Fabregat et al., 2017), and Biocarta (ver.2009_05_12) (Nishimura, 2001), and Gene Ontology annotations (Ashburner et al., 2000) were employed to represent associated biological processes. Statistical significance was defined by a p-value cut-off < 0.05 for all functional enrichment analyses. Each p-value was subsequently converted to a z-score by using the inverse normal cumulative distribution.

The prognostic analyses were performed using two datasets, GSE17536 (Smith et al., 2010) and TCGA-COADREAD containing patient survival data. The dataset of GSE17536 consists of 177 samples (24 patients with stage I, 57 patients with stage II, 57 patients with stage III, and 39 patients with stage IV) obtained from a patient cohort having the average age of 65.5 ± 13.1, whereas the dataset of TCGA COADREAD consists of 466 samples (89 patients with stage I, 170 patients with stage II, 130 patients with stage III, 62 patients with stage IV, and 15 patients with no stage information) acquired from a patient cohort having average age of 66.9 ± 12.5.

Survival analyses were performed by stratifying patients into high- and low-risk groups based on their prognostic index (PI), which is the linear component of the Cox model (PI = β₁x₁ + β₂x₂ + … + β_px_p, where β _p is the coefficient obtained from the Cox fit, x_p is the expression value of each gene) to investigate the prognostic performance of each DIP module. Analyses were performed with the web server SurvExpress (Aguirre-Gamboa et al., 2013) using datasets with clinical data. Signatures of survival in each risk group were estimated by Kaplan–Meier curves and Hazard Ratios (HR). The statistical significance of each curve was assessed by the cut-off for the log-rank p-value < 0.05. The hazard ratio (HR = (O1/E1)/(O2/E2)) was calculated to discover the significance of the survival curves based on the ratio between the relative death rate in group 1 (O₁/E₁) and the relative death rate in group 2 (O₂/E₂), where O denotes the observed number of deaths, and E denotes the expected number of deaths.

GeneXpharma (Turanli et al., 2017b), which is a publicly available platform presenting 50,304 gene-drug interactions among 4,344 genes and 11,939 drugs and employing statistical tests for the disease-gene-drug triad, was used for the network-based DR considering protein targets CDKN2A, GSK3B, HDAC2, and PML. A hypergeometric probability test was used to associate drugs with target proteins, and simulation results with p < 10⁻³ were accepted as statistically significant.

The 3-D structures of target proteins were obtained from Protein Data Bank (PDB) (Berman et al., 2002). The PDB identifiers were 1DC2 (CDKN2A) (Yuan et al., 2000), 6Y9R (GSK3B) (Buonfiglio et al., 2020), 6WBZ (HDAC2) (Yu et al., 2020), and 5YUF (PML) (Wang et al., 2018). Ligand binding sites of proteins were determined by PDBe (Velankar et al., 2016) and BIOVIA Discovery Studio (ver.20.1) (BIOVIA, 2016). The structures of candidate drugs were taken from PubChem (Kim et al., 2019). Autodock Vina software (ver.1.1.2) (Trott and Olson, 2019) was used for molecular docking analyses. All water molecules were deleted, and all polar hydrogens were added to the structure of a protein in the preparation of macromolecules. Molecular docking was performed five times for each ligand and each simulation yielded twenty poses. Exhaustiveness was set to eight for all docking calculations. Binding affinities were estimated to determine the importance of binding between protein targets and drug candidates. The top-scoring pose (with the lowest binding free energy) was selected for further analysis.

Abacavir (Selleckchem, S5215), Ribociclib (Selleckchem, S7440), Exemestane (Selleckchem, S1196), Voriconazole (Selleckchem, S1442), Bromocriptine mesylate (R&D systems, 0427/50), Tolcapone (Selleckchem, S4021) were dissolved in dimethyl sulfoxide (DMSO) and then diluted in 1x PBS. For each treatment step, 100 mM stock solution was diluted to a final concentration. Total DMSO concentration did not exceed 0.1% and control groups (CTRL) were also treated with the same concentration of DMSO. Nortriptyline hydrochloride (Selleckchem, S3698) and Theophylline (Sigma, T1633) were dissolved in 1x PBS.

CCD-841-Con cells (ATCC, CRL-1790, human healthy colon epithelial cell line) and HCT-116 cells (ATCC-CCL-247, human colorectal carcinoma cell line) were cultured in Dulbecco’s modified Eagle’s medium and Eagle’s Minimum Essential Medium, both supplemented with penicillin (100 units/ml), streptomycin (100 g/ml), and 10% fetal calf serum in a humidified atmosphere of 5% CO₂ and 95% air at 37°C. Later, cells were treated with varying abacavir, voriconazole, nortriptyline hydrochloride, tolcapone, and theophylline concentrations (10, 50, 100, 150, and 200 µM) in the first phase. Another wide exposure dose was used as 10, 50, 100, and 200 nM and 200–500 µM for bromocriptine mesylate, budesonide, exemestane, voriconazole, and theophylline. Also, nortriptyline hydrochloride concentration in the range of 10–50 µM was used to identify IC₅₀. Cell viability was determined using MTT (3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyl-2H-tetrazolium bromide) reduction by viable cells following the exposure durations of 24, 48, and 72 h 2 µl of MTT solution was mixed with 100 μL medium and incubated for 3 h at 37°C and 5% CO₂. At the end of this incubation period, the medium with MTT was discarded and 100 μL DMSO was added. The formazan crystals were dissolved by shaking at 150 rpm for 10 min. The intensity (OD) of the color formed was measured with a spectrophotometer with a microplate reader at a wavelength of 590 nm (Reference filter; 660 nm). Calculations were performed related to the absorbance of control samples which was equaled to 100%.

Two comprehensive transcriptome datasets associated with colorectal adenomas and tumors were recruited to apply differential interactome methodology for predicting high probability PPIs in tumor states and identifying differential PPIs. Two datasets were analyzed independently, and common dPPI signatures were considered in further analyzes. In the current formulation, upper and lower bounds (0.90 and 0.10, respectively) were used for q-values which represent the probability estimates. PPIs that had a normalized frequency of observation in either the normal or tumor phenotype greater than 20% were considered as significant. As a result, a total of 2,214 differential PPIs (dPPIs) were identified in the GEO dataset and 1,625 dPPIs were identified in the TCGA dataset. dPPIs were further classified as “significantly repressed in tumor phenotype” (if q < 0.10) or “significantly activated in tumor phenotype” (if q > 0.90). To this end, among the dPPIs, 718 interactions were repressed and 1707 interactions were activated in the GEO dataset, while 81 interactions were repressed and 1,557 interactions were activated in the TCGA dataset (Figure 1A).

The scale-free topology of the differential interactome network highlights the presence of hubs termed differentially interacting proteins (DIPs), indicating major changes in their interaction patterns during tumorigenesis (i.e., transition from “normal” to “tumor” phenotypes) (Ayyildiz et al., 2017). In this study, we identified DIPs and constructed a module around each DIP with its interaction partners, and the corresponding module was named adding the letter “m” to the beginning of the central DIPs (e.g., mMYC). The number of modules with more than 5 PPIs was 86 in TCGA and 177 in the GEO dataset. In both datasets, modules with at least 5 PPIs were considered and filtered based on their common interactions only. The number of common modules in two datasets (GEO and TCGA) were 24 (Figure 1B, Supplementary Table S1).

Among these modules, mMYC was the largest module with 116 members. Besides, the MYC protein was the most common protein observed in other modules. Moreover, mPSMB5, mCCT7, mRAN, mAPEX1, mHSD17B10, and mHDGF were also modules with more than 10 interactions. mPSMB5 had 33 interactions, as well as 30 interactions for the mCCT7, and mRAN had 29 interactions. For others, the number of interactions was 15 and 11, respectively.

Considering the members of the modules, gene set over-representation analyses were performed to identify biological processes, molecular pathways, and cancer hallmarks associated with the DIPs (Figure 2, Supplementary Table S1). The members of the modules were significantly enriched with several biological processes such as deubiquitination, protein modifications, cell cycle, neddylation, and transcription; signaling pathways such as PI3K-Akt, NOTCH, and Wnt signaling; cancer pathways such as transcriptional misregulation of cancer and microRNAs in cancer; and cancer pathways associated with oncogenic viruses such as Hepatitis B, Human papillomavirus (HPV), Epstein-Barr virus, and Human T-cell leukemia virus 1 (Figure 2A, Supplementary Tables S2, S3). Furthermore, the members of the modules were significantly associated with all cancer hallmarks (Senga and Grose, 2021), i.e., sustaining proliferative signaling, inducing angiogenesis, resisting cell death, deregulating cellular energetics, evading growth suppressors, activating invasion and metastasis, enabling replicative immortality, avoiding immune destruction, and genome instability and mutation (Figure 2B).

The evaluation of prognostic performances of DIP-centered modules was performed via Kaplan-Meier survival analyses through partitioning low- and high-risk groups regarding the expression levels of module genes. To cross-validate the results, the survival analyses were independently evaluated using two different datasets. As a result, prognostic performances of six modules (namely, mAPEX1, mCCT7, mHSD17B10, mMYC, mPSMB5, and mRAN) were found significant in both CRC datasets considering log-rank test p-values (p < 0.01) and hazard ratios (HR > 2) (Figure 3).

To identify potential drug targets, modules were further evaluated considering their sizes, members, prognostic performances, and enriched processes, pathways, and cancer hallmarks. Considering its relevance, we focused on the interactions of MYC in the module and filtered the proteins according to available 3-D structure information in PDB (Berman et al., 2002), and drug-gene interaction information in geneXpharma (Turanli et al., 2017b). The structures of 60 proteins out of 116 proteins were found in PDB, and drug interactions were identified for 39 proteins with available structural information. Durg repositioning simulation through geneXpharma resulted in 651 drugs associated with 39 target proteins. Then, these drugs were filtered according to FDA-approval, association with CRC in literature, and having the 3-D structures on PubChem (Kim et al., 2019). Finally, we identified eight drugs (abacavir, voriconazole, exemestane, nortriptyline hydrochloride, tolcapone, bromocriptine mesylate, ribociclib, and theophylline) targeting 4 proteins (PML, GSK3B, CDKN2A, HDAC2) (Figure 4).

In silico validation studies were carried out to evaluate the potential candidate drugs for the therapeutic strategies in CRC prior to in vitro viability assays. For this purpose, the 3-D structures of target proteins, i.e., 1DC2 (CDKN2A), 6Y9R (GSK3B), 6WBZ (HDAC2), and 5YUF (PML), and candidate drugs were retrieved from PDB (Berman et al., 2002) and PubChem (Kim et al., 2019), respectively, and the binding affinities of the drugs to their targets were estimated via protein-ligand molecular docking simulations using Autodock Vina (Trott and Olson, 2019) (Figure 5).

The binding affinity of the positive control case (theophylline) was predicted to be −5.1 kcal/mol. The most significant binding affinities were obtained for bromocriptine mesylate, nortriptyline hydrochloride, tolcapone, and exemestane, which were predicted as −9.7, −8.2, −7.3, and −7 kcal/mol, respectively. Similarly, the binding affinities of ribociclib, abacavir, and voriconazole were −6.9, −5.9, −5.7 kcal/mol, respectively. Besides, the binding affinity of nortriptyline hydrochloride which targets PML protein was -6.1 kcal/mol. According to the molecular docking simulations, all potential drugs (bromocriptine mesylate, tolcapone, nortriptyline hydrochloride, ribociclib, exemestane, voriconazole, and abacavir) display significantly higher binding affinities when docked to their target proteins, compared to the positive control case (Figure 5). The details of the top-scoring pose of five docking simulations for each ligand are displayed in Supplementary Figure S1.

To have an insight into the antitumor potential of drugs for CRC treatment, in vitro cell viability assay was performed in HCT-116 (carcinoma) cell line. Concentration ranges as 10–200 uM for abacavir and tolcapone, 5–100 uM for ribociclib, 200–500 uM for exemestane, voriconazole, and bromocriptine mesylate, 10–50 µM for nortriptyline hydrochloride, and 20–500 µM for theophylline were tested to find out the IC₅₀ values in HCT-116 cells.

After 24, 48, and 72 h drug treatments, IC₅₀ values were determined as 100 μM at 72 h for abacavir, 100 μM at 72 h for ribociclib, 500 μM at 48 h for exemestane, 40 μM at 48 h for nortriptyline hydrochloride, 500 μM at 72 h for theophylline and 10 μM at 72 h for tolcapone (Figure 6). Bromocriptine mesylate and voriconazole did not decrease viability at tested concentrations and time points.

In addition, the toxicity of the drugs was tested in and CCD-841-Con healthy epithelial cells in the obtained IC₅₀ values. Ribociclib was the only drug that caused death in healthy cells when compared to control (Figure 7).