Next Article in Journal
PCAT6 May Be a Whistler and Checkpoint Target for Precision Therapy in Human Cancers
Next Article in Special Issue
RKIP Pleiotropic Activities in Cancer and Inflammatory Diseases: Role in Immunity
Previous Article in Journal
Childhood Malignant Brain Tumors: Balancing the Bench and Bedside
Previous Article in Special Issue
RKIP: A Pivotal Gene Product in the Pathogenesis of Cancer
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Cancers
Volume 13
Issue 23
10.3390/cancers13236098
cancers-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Functional Network Model of the Metastasis Suppressor PEBP1/RKIP and Its Regulators in Breast Cancer Cells

by
Mahmoud Ahmed
†,
Trang Huyen Lai
†,
Wanil Kim
and
Deok Ryong Kim
*
Department of Biochemistry and Convergence Medical Science, Institute of Health Sciences, Gyeongsang National University College of Medicine, Jinju 527-27, Korea
*
Author to whom correspondence should be addressed.
†
These authors contributed equally to this work.
Cancers 2021, 13(23), 6098; https://doi.org/10.3390/cancers13236098
Submission received: 25 October 2021 / Revised: 26 November 2021 / Accepted: 30 November 2021 / Published: 3 December 2021
(This article belongs to the Special Issue RKIP: A Pivotal Gene Product in the Pathogenesis of Cancer)
Browse Figures
Review Reports Versions Notes

Abstract

:

Simple Summary

We extracted known interactions between metastasis suppressors and their regulators from scientific studies. Next, we used publically available data of drug treatments to identify which of them potentially perturbs these interactions. Finally, we studied the effect of several of these drugs on a particular metastasis suppressor called RKIP and found that our model accurately predicted its regulation in breast cancer cells. Our approach can discover alternative mechanisms of existing cancer drugs and repurpose them in different disease types.

Abstract

Drug screening strategies focus on quantifying the phenotypic effects of different compounds on biological systems. High-throughput technologies have the potential to understand further the mechanisms by which these drugs produce the desired outcome. Reverse causal reasoning integrates existing biological knowledge and measurements of gene and protein abundances to infer their function. This approach can be employed to appraise the existing biological knowledge and data to prioritize targets for cancer therapies. We applied text mining and a manual literature search to extract known interactions between several metastasis suppressors and their regulators. We then identified the relevant interactions in the breast cancer cell line MCF7 using a knockdown dataset. We finally adopted a reverse causal reasoning approach to evaluate and prioritize pathways that are most consistent and responsive to drugs that inhibit cell growth. We evaluated this model in terms of agreement with the observations under treatment of several drugs that produced growth inhibition of cancer cell lines. In particular, we suggested that the metastasis suppressor PEBP1/RKIP is on the receiving end of two significant regulatory mechanisms. One involves RELA (transcription factor p65) and SNAI1, which were previously reported to inhibit PEBP1. The other involves the estrogen receptor (ESR1), which induces PEBP1 through the kinase NME1. Our model was derived in the specific context of breast cancer, but the observed responses to drug treatments were consistent in other cell lines. We further validated some of the predicted regulatory links in the breast cancer cell line MCF7 experimentally and highlighted the points of uncertainty in our model. To summarize, our model was consistent with the observed changes in activity with drug perturbations. In particular, two pathways, including PEBP1, were highly responsive and would be likely targets for intervention.

1. Introduction

Several approaches for screening and identifying bioactive compounds exist and are being used to find effective remedies for cancer. Most strategies focus on quantifying the phenotypic effects of different biological models [1]. High-throughput technologies like massive parallel sequencing and mass spectrometry have the potential to augment these strategies by explaining the effect of panels of drugs on biological entities such as gene expression and protein abundance. Reverse causal reasoning enables the use of quantitative measurements of entities in a given pathway to make claims about their functions [2,3]. This method integrates existing biological knowledge and data of gene and protein abundance to infer their function [4]. In addition to taking advantage of prior knowledge, this approach utilizes the network structure, which is key in biological pathways and avoids making assumptions about equating gene expression with protein abundance or activity [5,6].
Breast cancer is the most prevalent type of cancer among women, and the vast majority of deaths are due to metastasis [7,8]. Cancer metastasis is a complex process that involves disseminating cancer cells from primary solid tumors to distant parts of the body, where they form secondary tumors. Effective cancer therapies require an understanding of the underlying mechanisms and regulatory pathways of genes involved in metastasis. A group of genes, known as metastasis suppressor genes (MSG), prevents the formation of overt metastases [9]. Although multiple MSGs have been discovered, the detailed mechanism by which they overcome metastasis is not well understood. One of the metastasis suppressors is phosphatidylethanolamine binding protein (PEBP1), also known as RAF kinase inhibitory protein (RKIP) [10]. PEBP1 inhibits cancer metastasis at different stages, including epithelial to mesenchymal transition (EMT), migration, and invasion [11,12,13]. Low expression of PEBP1 is associated with cancer metastasis and poor prognosis [14].
Here, we attempted to integrate two sources of data in the form of biological knowledge and gene expression profiles to build a functional model of metastasis in breast cancer. We used a reverse causal reasoning approach to infer the function of crucial metastatic proteins from the perturbation of expression of other gene products known to interact with them. We evaluated the model for agreement in the direction of change with drug treatments that inhibit the growth of the cancer cells. Finally, we prioritized a few key pathways involving the anti-metastasis PEBP1 for further testing.

2. Materials and Methods

2.1. Protein-Protein Interactions, Literature Search, and Biological Expression Language

We used the String database to identify protein–protein interactions (PPI) among a list of seventeen MSGs (Table 1), and eleven transcription factors (TFs) (Table 2) [15]. Most of the retrieved interactions were based on text mining, where two or more entries co-occur in published articles. We surveyed the linked studies for evidence to establish causal relations between the entities. We encoded these interactions in the biological expression language (BEL) and assigned at least one reference to each (Table 3) [16]. In BEL, two entities (e.g., r, mRNA or p, protein) interact when a subject produces an effect (increase, decrease, directlyIncrease or directlyDecrease) on the object. We used these interactions to build the functional layer of the metastatic suppressors and their regulators.

2.2. Knockdown of Metastasis Suppressors and Transcription Factors in MCF7

We used two datasets where each MSG and TF was knocked down in the MCF7 breast cancer cell line. The first is a dataset that consists of 77 microarray samples [28]. Seventeen MSGs were knocked down (siRNA) separately in MCF7 (Table 1). The control set consisted of wild-type (0 h) and cells treated with scramble siRNA at 48 or 96 h. We obtained the processed data using GEOquery [29]. We mapped the probe intensities to gene symbols and applied differential expression between knockdown and control using limma [30]. We compiled a second dataset from the KnockTF database for the eleven TFs knockdown in the same cell line (Table 2) [31]. KnockTF curates gene expression profiles of TFs knockdown from the gene expression omnibus (GEO) and presents the results in the form of differential expression output (fold-change and p-values).

2.3. Pharmacological Perturbations of Breast Cancer Cell Lines

Two sources of publicly available drug perturbation data were used in this study. The first is growth inhibition (GRmax) data of drug perturbations (n = 35) in breast cancer cell lines (n = 71, including MCF7) [32]. Adenosine triphosphate content was measured as a proxy for the cell count in drug treated cell cultures relative to DSMO-treated controls. We included in the analysis only the drugs with at least one effective dose. This dataset was used as evidence for the efficacy of these drugs in producing phenotypic effects on cancer cells. The second data source is the library of integrated network-based cellular signatures (LINCS), which contains gene expression profiles of multiple cells under different types of perturbations, including compounds, gene overexpression, knockdown, or knockouts [33]. Several cell lines were profiled for gene expression using the L1000 technology, which measures the abundance of 1000 landmark genes and infers the expression of the rest of the genes. We used a subset of drug perturbations (n = 35, same set of drugs as above) in cancer cell lines (n = 67) that were also surveyed for drug sensitivity. We obtained the expression profiles of all genes (level 3 of the LINCS data) using Slinky [34].

2.4. Network Perturbation Amplitude (NPA)

Network perturbation amplitude (NPA) models the perturbation of a causal network as the changes in the expression of connected downstream nodes [4]. This method assumes that the downstream gene expression reflects the function of the entities in a biological pathway (e.g., the activity of a given TF). The input to this method is a two-layer network model of the pathway of interest and expression profiles of two or more conditions. First, the functional layer (backbone) encodes the causal relations between the entities of interest. The second is a transcription layer (terminal) with all the nodes downstream from each node in the backbone. We used the gene expression profiles from LINCS pharmacological perturbations to score the model. NPA is an R package that efficiently implements the method with the same name [35].
Table
Table 3. Interactions of metastasis suppressor genes and their transcription factors.
Table 3. Interactions of metastasis suppressor genes and their transcription factors.
SubjectObjectRef.Interaction
Interactions between metastasis suppressors
CD44CASP8[36]p(CD44) decreases act(p(CASP8))
CD82CD44[37]p(CD82) decreases p(CD44)
CD82BRMS1[38]p(CD82) increases r(BRMS1)
CD82CDH1[39]p(CD82) increases r(CDH1)
CDH1CASP8[40]p(CDH1) increases act(p(CASP8))
GSNCDH2[41]p(GSN) increases r(CDH2)
MAP2K3MAPK14[42]act(comp(p(MAP2K6), p(MAP2K3)), ma(kin)) directlyIncreases p(MAPK14)
MAP2K4MAP2K7[43]p(MAP2K4) decreases act(p(MAP2K7), ma(kin))
MAP2K4MAPK14[42]p(MAP2K4) directlyIncreases act(p(:MAPK14), ma(kin))
MAP2K4CASP8[44]p(MAPK14) increases p(CASP8)
MAP2K6MAPK14[42]p(MAP2K6) directlyIncreases p(MAPK14)
MAP2K7CDH2[45]p(MAP2K7) increases r(CDH2)
MAPK14CDH1[46]p(MAPK14) increases p(CDH1)
MAPK14CASP8[47]act(p(MAPK14), ma(kin)) directlyDecreases p(CASP8)
MAPK14SAP1A[48]act(p(MAPK14), ma(kin)) directlyIncreases p(ets-Domain Protein Elk-4)
MAPK14RUNX2[49]p(MAPK14) increases r(RUNX2)
NME1BRMS1[50]p(NME1) decreases r(BRMS1)
NME1AKAP12[50]p(NME1) increases r(AKAP12)
NME1PEBP1[51]r(NME1) increases r(PEBP1)
PEBP1MAP2K6[52]p(PEBP1) increases act(p(MAP2K6))
PEBP1MAP2K3[52]p(PEBP1) increases act(p(MAP2K3))
PEBP1MAPK14[52]p(PEBP1) increases act(p(MAPK14))
PEBP1CDH1[53]p(PEBP1) increases r(CDH1)
RUNX2CDH2[54]p(RUNX2) decreases r(CDH2)
RUNX2CDH1[54]p(RUNX2) increases r(CDH1)
TGFB1CDH1[55]p(TGFB1) decreases r(CDH1)
TGFB1CDH2[56]p(TGFB1) decreases r(CDH2)
Interactions between metastasis suppressors and regulators
CDH1YBX1[57]p(CDH1) decreases p(YBX1)
CDH1FOXM1[58]p(CDH1) decreases p(FOXM1)
ESR1CDH2[59]p(ESR1) decreases r(CDH2)
ESR1CDH1[60]p(ESR1) directlyIncreases r(CDH1)
ESR1BRMS1[61]p(HGNC:ESR1) directlyIncreases r(HGNC:BRMS1)
ESR1MTA3[62]p(ESR1) increases r(MTA3)
ESR1RUNX2[63]p(ESR1) increases r(RUNX2)
FOSCDH1[55]p(FOS) decreases r(CDH1)
FOSRUNX2[64]p(FOS) decreases r(RUNX2)
FOSTGFB1[55]p(FOS) directlyDecreases r(TGFB1)
FOSCD44[65]p(FOS) increases r(CD44)
GATA3CD44[66]r(GATA3) directlyIncreases r(CD44)
GATA3CD44[66]r(GATA3) directlyIncreases r(CD44)
GATA3CD44[66]p(GATA3) increases r(CD44)
GATA3CD44[66]p(GATA3) increases r(CD44)
HIF1ACD44[67]p(HIF1A) increases r(CD44)
MAP2K4FOS[68]p(MAP2K4) increases p(FOS)
MAPK14FOS[48]act(p(MAPK14), ma(kin)) increases p(FOS)
MAPK14GATA3[69]p(MAPK14) increases act(p(GATA3))
MTA3CDH2[62]p(MTA3) decreases r(CDH2)
MTA3CDH1[62]p(MTA3) increases r(CDH1)
RELABRMS1[70]p(RELA) directlyDecreases r(BRMS1)
RUNX2HIF1A[71]p(RUNX2) increases act(p(HIF1A))
RUNX2SPDEF[72]p(RUNX2) directlyDecreases r(SPDEF)
SAP1AFOS[48]p(ets-Domain Protein Elk-4) directlyIncreases r(FOS)
SATB1NME1[73]p(SATB1) directlyDecreases r(NME1)
SATB1BRMS1[73]p(SATB1) directlyDecreases r(BRMS1)
SATB1CD82[73]p(SATB1) directlyDecreases r(CD82)
SATB1CDH1[73]p(SATB1) directlyDecreases r(CDH1)
SNAI1CDH1[74]p(SNAI1) directlyDecreases r(CDH1)
SNAI1PEBP1[75]p(SNAI1) directlyDecreases r(PEBP1)
SNAI1CASP8[76]p(SNAI1) decreases act(p(CASP8))
TFAP2CCD44[77]p(TFAP2C) directlyDecreases r(CD44)
YBX1CD44[78]p(YBX1) directlyIncreases r(CD44)
Interactions between transcription factors
ESR1HIF1A[79]p(ESR1) directlyIncreases r(HIF1A)
ESR1GATA3[20]p(GATA3) increases p(ESR1)
FOXM1GATA3[80]p(FOXM1) directlyDecreases r(GATA3)
MTA3SNAI1[81]p(MTA3) directlyDecreases r(SNAI1)
NR5A2ESR1[82]p(NR5A2) directlyIncreases r(ESR1)
RARAFOS[83]p(RARA) decreases act(p(FOS))
RELASNAI1[84]p(RELA) directlyIncreases r(SNAI1)
TFAP2CESR1[85]p(TFAP2C) directlyIncreases r(ESR1)
TNFSF10CASP8[86]p(TNFSF10) increases p(CASP8)

2.5. Measures of Agreement

We defined two notions of agreement between expectations and observations of the effects of drug perturbations on the metastasis network model. For every node (u), we created a graph of all n nodes connected to it by an edge (e) (subnetwork). Given a binarized perturbation coefficient of drug treatment, we formed an expectation by multiplying +1 or −1 with the direction of interaction +1 or −1 for every edge, u × e . We compared the expected sign to the observed perturbation coefficient of the downstream nodes ( x ′ ). The average for every subnetwork we termed a concordance rate. We similarly computed an agreement estimate for all possible directed graphs that lead to a particular node (path). The difference from above is that the upstream node is not fixed, but changes as we move one edge down the path. The average agreement between the observed and expected effect for every path we termed a coherence rate. The following formula shows how the agreement rates were calculated.
1 n ∑ i = 1 n x i x = 1 , if u . e = x ′ 0 , otherwise u = 1 , if u is activated − 1 , if u is repressed
where x is a component connecting a node u to a downstream node by an edge e and x ′ is the observed effect of drug perturbation on that downstream node in the component of the subnetworks (concordance) or paths (coherence).
In both cases, the agreement was expressed in terms of Cohen’s κ [87]. This is a value between −1 (worse) and 1 (better) while taking into account the expected agreement by chance. The following shows how to calculate it.
κ = ( observed agreement − expected agreement ) / ( 1 − expected agreement )
Finally, we compared the probability distribution of the agreement measures to randomly generated values using the Kolmogorov–Smirnov (KS) test. D + is the maximum distance between the cumulative distribution functions (ECDF).

2.6. Reagents and Drugs

Reagents and drugs utilized in this study were purchased as follows: RPMI-1640 media (11875-119), fetal bovine serum (FBS; 16000-044), Epirubicin hydrochloride (CAS 56390-09-1), Vorinostat/SAHA (CAS 149647-78-9), Methotrexate hydrate (CAS133073-73-1), Cisplatin (CAS 15663-27-1), Sorafenib (CAS 284461-73-0) from Sigma-Aldrich (St. Louis, MO, USA); Imatinib (CAS 220127-57-1) from STEMCELL Technologies Inc. (Vancouver, BC, Canada); Trizol reagent (15596026) from Invitrogen (Carlsbad, CA, USA); DNase I Solution (1 unit/ μ L), RNase-free (89836), and RevertAid First Strand cDNA Synthesis Kit from Thermo Scientific (Waltham, MA, USA); amfiSure qGreen Q-PCR Master Mix(2X), Without ROX (Q5600-005) from GenDEPOT (Katy, TX, USA).

2.7. Cell Culture and Drugs Treatment

MCF-7 breast cancer cells were cultured in RPMI-1640 supplemented with 10% Fetal Bovine Serum (GIBCO) and 100 μ g/mL streptomycin and incubated in a 37 ∘ C humidified atmosphere containing 5% CO 2 . For drug treatment, cells ( 3 × 10 5 cells) were plated in a 6-well plate and further incubated for 24 h. Then, cells were treated with indicated drugs: Epirubicin, Vorinostat, Methotrexate, Cisplatin, Sorafenib, and Imatinib at 1.5 μ M concentration and collected for assays after 24 h of treatment.

2.8. RNA Extraction, and RT-qPCR

Total RNAs were extracted from MCF7 breast cancer cells using Trizol reagent according to the manufacture’s manual. Prior to cDNA synthesis, RNA samples were treated with DNase I solution to remove trace amounts of DNA. The prepared RNA samples were then used as templates for reverse transcriptase to generate First-strand cDNA using Thermo Scientific RevertAid First Strand cDNA Synthesis Kit. The first-strand cDNA synthesis products were used directly in qPCR using the amfiSure qGreen Q-PCR Master Mix kit. Primers of target genes used in the assay are presented in Table 4. The expression of the five genes was quantified in the treated samples relative to the control gene GAPDH and the control condition of the DMSO treatment. Δ Δ C t model was applied using the pcr R package [88]. Student t-test was used to compare the relative expression in each treatment to the control DMSO treated cells. p-values < 0.05 were considered significant. Experiments were performed in five or more replicates.

2.9. Software Environment and Reproducibility

This analysis was performed in R and using Bioconductor packages [89,90]. The software environment was packaged and distributed as a Docker image (https://hub.docker.com/r/bcmslab/antimetastatic, accessed on 29 November 2021). The code to run the analysis and reproduce the figures and tables in this manuscript is available as open-source (GPL-3) (https://github.com/BCMSLab/antimetastatic, accessed on 29 November 2021).

3. Results

3.1. A Workflow for Building a Network of the Metastasis Suppressors and Their Regulators

In this study, we built and evaluated a network model of MSGs and their transcriptional regulators in breast cancer (Figure 1A). First, we identified PPIs of seventeen MSG and eleven TFs in the String database. Then, we vetted and supplemented every link with information from the relevant literature to determine the type of entity, direction of interaction, and evidence. We added a few interactions that were not present in the PPIs.
Since these causal links originated from experiments in different organisms and conditions, we further filtered the model to keep only the interactions relevant to the context of breast cancer. We obtained and analyzed the knockdown datasets of each of the input genes and regulators in the MCF7 cell line. We resolved the conflicts between the reported interactions and fold-change from the knockdown data by prioritizing the context of the interactions, the strength of the evidence, and the effect size (Figure 1B). We used differential expression to identify possible links between the nodes that were not previously reported and add them to the network. We compiled the results were in a context-specific metastasis model using BEL. The full list of included regulatory links is reported in Table 3 into three categories: within MSGs, within TFs, and in between the two groups of gene products.
Finally, we evaluated the context-specific model using gene expression and growth rate data in MCF7 treated with different compounds. The context-specific interactions served as a backbone to a transcript layer derived from the knockdown data to form a two-layer network. We scored the model using NPA under different conditions. We then compared the effect of the drugs on the entire network to their impact on the growth rate. Next, and given the expected effect of drugs on each node, we compared this predicted change on the downstream nodes to the observed using two different measures of agreement (Figure 1C,D). Moreover, we constructed a reliable model of highly consistent paths for further exploration.

3.2. Identifying Possible Interactions of Metastasis Suppressors Using Text Mining of the Literature

We took as a starting point the known protein–protein interactions between seventeen MSGs (Table 1) and their eleven TFs (Table 2). The initial query on the String database resulted in 93 interactions. Most interactions were based on text mining or the proximity of terms (gene products) in the text of one or more published study. We manually examined the provided evidence for each interaction. We also filtered, and supplemented with information from related studies. The interactions were coded in BEL and were used as the functional layer for the NPA analysis (Table 3). Each pair of interacting gene products were considered as a subject and an object. The type and direction of the relation between the two entities was represented according to the BEL standards.

3.3. Contextualizing Metastasis Suppressor Interactions in Breast Cancer Cells

Inevitably, conflicts arise between the relations suggested in the literature and the manifest changes of knocking down some genes on the expression of others. For example, not all interactions are relevant in the biology of breast cancer because they come from studies of different conditions. Additionally, not all interactions are reported in previous studies. We used a data-driven approach to specify the associations that are most relevant and augment the network with interactions based on changes in gene expression. The efficiency of the knockdown was verified by checking the expression of the corresponding genes (Figure 2, bottom). Most knocked down genes had a significant fold-change (log 2 FC < −1 and p-value < 0.05) when compared to the control cells. Similarly, a wide-ranging change (between −4 and 2 log 2 FC) in expression in the metastasis suppressors and regulators resulted from the knockdown, in most cases (Figure 2, top).
We compared the observed fold-change to the directions of previously curated interactions. We kept the ones that had no evidence to the contrary as a result of the knockdown of the subject. We removed the links with a significant (absolute log 2 FC > 1 and p-value < 0.05) change in the opposite direction of change when the upstream node was knocked down. Finally, we coded the significant effects of the knockdown of one node (subject) on others (object) as an interaction between the mRNA to directly increase or decrease. The resulting network contained 78 edges between the input 35 unique nodes (Table 3). The recorded interactions included relations between the coding RNA, protein, and protein activities. That direction of interaction was either positive (increase/directlyIncrease) or negative (decrease/directlyDecrease) (Figure 3A).
TFs seem to have extensive interactions, reflecting the hierarchical nature of the regulator networks reported in other studies. Only a few nodes in the network had many edges, and most nodes had, on average, only two edges (Figure 3B, top). The highly connected nodes were either TFs that received signals from other regulators, key nodes in signaling pathways, or effector nodes that are tightly regulated. Examples include FOS and ESR1 as TFs, MAP Kinase subunits as signaling nodes, or CDH1 as an effector. In addition to filtering in the relevant interactions, we used the knockdown dataset to build a transcript layer for the functional model. The numbers of nodes connected to those in the functional layer were more normally distributed. Similarly, the distribution of up and down-regulation is normal but down-regulated nodes were more frequent (Figure 3B, bottom).

3.4. Evaluating the Metastasis Model Using Drug Perturbation Data

To evaluate the relevance of the metastasis network, we used a dataset of growth inhibition rates and gene expression profiles under drug perturbations. In this dataset, the growth rate was quantified in breast cancer cell lines (n = 71), including MCF7. Cell cultures were treated with several compounds (n = 35) and compared to the baseline replication rate. We included in the analysis only the drugs with at least one effective dose in terms of significantly (GRmax < 1) inhibiting cell growth (Figure 4A). Drug responses in most cell lines were correlated with the response in MCF7 (Figure 4B). NPA infers the activity of the nodes in the functional layer of the metastasis network using a transcript layer derived from the knockdown dataset. This method returns a perturbation coefficient that indicates the magnitude and direction of changes (−1, inhibition, and 1, activation) for the network as a whole and for each node in the graph (Figure 4C). The perturbation coefficients ranged from −0.1 to 0.2, with most nodes being significantly (confidence intervals not including zero) perturbed at least once. The perturbation amplitudes of 67 cancer cell lines correlated well with the NPA value in MCF7 (Figure 4D).
Next, we divided the metastasis graph into smaller modules and evaluated the agreement between the expected direction of change and the observed coefficients for each perturbation. We considered every node as an upstream connected to one or more nodes by one edge (subnetwork). The effect of the drugs on the upstream node was binarized (activation > 0, repression < 0) and used to infer the expected impact on the downstream nodes. We compared the expectations to the observed perturbation coefficients of the downstream node after similar binarization (concordance; see Section 2).
Higher concordance supports the consistency of the suggested model. Cohen’s κ was employed to evaluate the agreement between expectations and observation ( κ > 0, better and κ < 0, worse), taking into account the chance agreement (50%). The average concordance of each drug was higher (D + = 2.9 and p-value < 0.0001 in KS test) than expected by chance alone (Figure 5A). The network accurately predicted the direction of drug perturbations. Averaging for each node, a high concordance (D + = 2.6 and p-value < 0.0001) was also observed (Figure 5B). The correlation between the concordance of the subgraph and the effect of the drug on its source node, a potential source of bias, was very weak (r < −0.02 and p-value > 0.9). We observed higher than chance (D + = 2.47 and p-value < 0.0001) of consistency in other cell lines treated with the same set of drugs (Figure 5C).

3.5. Constructing a Model of PEBP1 and Its Interaction with Other Metastasis Suppressors

To demonstrate the utility of this model, we used the same perturbation data to prioritize key nodes and consistent edges between them. We isolated all paths in the graph leading to PEPB1 (n = 11). Paths are all interactions in the directed graph that end with the node corresponding to a protein. We multiplied the binarized perturbation coefficients of the upstream nodes by the sign of the edges between the nodes in every path. We compared these expectations were the observed perturbation coefficients on the downstream nodes of the paths (coherence; see Section 2). Not all drugs produced coherent effects on PEBP1 pathways, but the agreement was bigger than to be by chance (D + = 2.44 and p-value < 0.0001 in KS test) (Figure 6A). Most paths were coherent (D + = 1.45 and p-value < 0.0001) indicating the strong relevance of these pathways/interactions (Figure 6B). The main sources of bias (upstream coefficients) were not strongly relevant (r < −0.1 and p-value > 0.7). The observed directions of interaction were also consistent (D + = 3.01 and p-value < 0.0001) with predictions across different cancer cell lines when treated with the same drugs (Figure 6C).
We used the top five coherent (coherence > 0.6) PEBP1 paths to highlight interactions that are strongly relevant to breast cancer cell metastasis and are highly responsive to treatment with different drugs (Figure 7A). We constructed a model of regulatory interactions between PEBP1 and its upstream regulators with the top five coherent paths represented by dash-lines (Figure 7B). Interestingly, we found that the pathways flowed in two main directions to enhance or suppress PEBP1 activity—the upstream proteins of PEBP1 function as TFs/regulators. One involves RELA (p65) and SNAI1, which were previously reported to inhibit PEBP1 transcription. The other involves the estrogen receptor (ESR1), which induces PEBP1 through the kinase NME1, which to our knowledge has not been reported before.

3.6. Validating the Model Predictions

To test the validity of the predictive model of PEBP1/RKIP regulation, we selected several drugs that activate or repress PEBP1 and tested their effect on gene expression in MCF7. Indeed, four out of 6 drugs conformed to the expected direction of regulation. Epirubicin induced the expression of PEBP1 while Sorafenib, Cisplatin and Imatinib repressed it (Figure 8A). One repressor (Sorafenib) seems to produce its effect on PEBP1 through either of the two suggested pathways (Figure 8B). Sorafenib treatment induced RELA and SNAI1, which lowered PEBP1 and repressed ESR1 and NME1, which activated it. The effect size on NME1 was small and not significant at a cuttoff of p = 0.05 . By contrast, Cisplatin and Imatinib repressed PEBP1 through one or the other regulatory pathways. Cisplatin induced RELA and SNAI1 which inhibited PEBP1 (Figure 8D), while Imatinib repressed PEPB1 through ESR1 and NME1 inhibition (Figure 8D). The activator, Epirubicin, induced the expression of PEBP1 through activating NME1 (Figure 8E). Our model predicted the direction of regulation of PEBP1 and the pathway by which this regulation was achieved (Figure 8F).

4. Discussion

We used text mining datasets and a manual literature search to extract evidence for possible interactions between several metastasis suppressors and their regulators. We filtered these interactions in the context of breast cancer using a knockdown dataset in the cell line MCF7. The resulting interactions were coded in BEL to build a functional metastasis model. Finally, we used a reverse causal reasoning approach to evaluate and prioritize these interactions and extract pathways that are most consistent with drug treatments that inhibit cell growth.
High-throughput technologies such as massive parallel sequencing and mass spectrometry produce simultaneous measurements on many biological entities. Interpreting these measurements requires expert biological knowledge. An approach to formalize the use of knowledge and data is to use causal reverse reasoning [4]. Here, we employ such an approach to infer the function of metastatic suppressors using curated causal links and the perturbation data. The implicit advantage of this approach is also to make claims about the role of gene products from measurements of gene expression. Castro and colleagues used a similar approach to construct a model of inflammatory bowel disease (IBD). They were able show that three pathways that contribute to the disease were regulated differently in the two subtypes of IBD [91]. Another study identified a mechanism the mediates the difference between the subtypes in response to environmental exposure such as cigarette smoke [92].
Studies based on text mining or manual curation alone can suffer from several drawbacks. Some of the identified interactions can only be relevant in a particular biological or pathological context. In addition, inconsistencies arise with changing the context of the interactions. Finally, inevitable biases exist in the literature toward studying specific pathways or models. We attempted to address these issues in our study with a number of strategies, which we will discuss next.
We took all possible interactions as defined by text mining of relevant studies and manual curation as a prior (Figure 1A). Next, we evaluated the relevance and consistency of these links. Interactions that were not contradicted by the knockdown of their corresponding nodes are likely to be present in this context (Figure 1B). This filtering step increases the likelihood of the reported links being relevant. Similarly, consistency in the direction of change with drug perturbations increases our confidence in the relevance of some interactions which we initially included. We implemented intuitive consistency measures for all possible subnetworks (Figure 1C) and for selected paths (Figure 1D) that include PEBP1 as the target node.
In addition to filtering out irrelevant and insignificant interactions, we used the knockdown dataset to identify regulatory links that were not previously reported in the literature. We inferred potential edges from significant changes in gene expression as a result of the knockdown of a given node (Figure 2). These links would augment the model by increasing the number of edges and mitigating biases arising from tendencies in the literature toward studying certain pathways and specific disease contexts. We suspected that larger effect sizes (perturbation amplitudes) could bias the agreement measures. However, neither in the case of subnetwork concordance nor the path coherence was this bias big enough to influence the results (Figure 5C and Figure 6C).
The value of any predictive model is to make testable hypotheses. Here, we suggested two key pathways regulating PEBP1 (Figure 7). First, the inhibition of the zinc transcriptional repressor (SNAI1) results in PEBP1 suppression [75]. In contrast to PEBP1’s role in metastasis, SNAI1 induced breast cancer EMT and metastasis by directly repressing the epithelial markers E-cadherin (CDH1) [74]. Martin and colleagues found SNAI1 to be overexpressed in invasive metastatic breast cancer compared to normal breast tissue [93]. Moreover, RELA (P65), which is one of five subunits of the Nuclear factor- k a p p a B family, positively regulated SNAI1 transcription. NF- κ B binds directly to the SNAI1 promoter, resulting in increased SNAI1 transcription [94]. Overall, NF- κ B, SNAI1, and PEBP1 form a feedback loop that regulates the development of metastasis and resistance to apoptosis [95]. MTA3 is another upstream node of SNAI1 and suppresses its transcription. Similarly, MTA3 represses the transcription of a series of EMT-promoting genes such as ZEB2 and N-cadherin [81]. Our model is consistent with the above-mentioned studies.
Another pathway to regulate PEBP1 in our model is through ESR1 and NME1. NME1 dramatically inhibits the NF- κ B signaling pathway and alters the transcription of metastasis-related genes, reducing metastatic pulmonary cells [96]. Additionally, NME1 regulates gene expression in breast cancer cells. Moreover, NME1-dependent genes have a prognostic value as they predicted survival in breast cancer patients [50]. Here, we propose for the first time an alternative link between the ESR1 and PEBP1 induction via the kinase NME1.
We experimentally tested the effect of several drugs on PEBP1/RKIP and its regulators. In addition to verifying the predicted interactions, we wanted to highlight the uncertainty associated with this model. The model had three points of uncertainty. The first is whether a given drug produces the predicted effect on the terminal node of interest. Four out of the six tested drugs conformed to the expectations. The second uncertainty is whether that effect is achieved through both, either or a third unspecified pathway. Sorafenib inhibited PEBP1 by activating the inhibitory pathway of RELA/SNAI1 and inhibiting the other activation pathway of ESR1/NME. By contrast, Cisplatin and Imatinib repressed PEBP1 through only one of the two pathways. The third uncertainty stems from the fact that the farther back in the pathway we travel, the less accurate the predictions are. For example, Epirubicin activated PEBP1 by inducing NME1, but the effect on ESR1 was not as expected.
The validity of this model is limited to how much the included interactions represent the underlying biology. Another limitation is the amount of gene expression signal captured by the knockdown and perturbation data. We carefully curated existing studies to have as many and as accurate as possible links between metastasis suppressors. This model is guaranteed to evolve with the rapidly accumulating biomedical literature. As there is no existing high-throughput method to measure the activity of diverse biological entities in parallel, we relied on a reverse engineering approach to infer the activity from gene expression changes in their downstream nodes. The reverse reasoning approach assumes that the changes in gene expression of a group of nodes reflect the activation or inhibition of their upstream proteins.
The regulatory links in the final model were derived from the literature and prioritized by our analysis. However, further validation of these interactions is necessary first to test their relevance and the reliability of this approach. We supported our observations by devising two notions of consistency, concordance, and coherence. The first of the two shows the general agreement between expectations and observations in the subgraphs and the response to drug treatment. On the other hand, coherence shows consistency in particular pathways and can be thought of as a measure of the relevance of specific links. We elected to use those two measures of agreement for ease of interpretation rather than more sophisticated notions of consistency. We carried out experimental validation and found that our model is most reliable in predicting the drug effect on a target node, but uncertainty increases as to through which pathway and how many upstream nodes are regulated.

5. Conclusions

In conclusion, we suggested that the metastasis suppressor PEBP1 is on the receiving end of two key regulatory mechanisms. One inhibitory pathway involves RELA (NFKB p65 subunit) and SNAI1, previously reported to interact with PEBP1. The other pathway involves the estrogen receptor (ESR1), which induces PEBP1 through the kinase NME1. Those two pathways were highly responsive to pharmacological perturbations and would be likely targets for intervention against metastatic breast cancer.

Author Contributions

M.A., conceptualization, methodology, software, formal analysis, visualization, and writing. T.H.L., conceptualization, methodology, data curation, data generation and writing. W.K., methodology, writing. D.R.K., conceptualization, writing, supervision, project administration, and funding acquisition. All authors have read and agreed to the published version of the manuscript.

Funding

This study was supported by the National Research Foundation of Korea (NRF) grant funded by the Ministry of Science and ICT (MSIT) of the Korean government [2015R1A5A2008833 and 2020R1A2C2011416].

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data are contained within the article.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
MSGMetastasis suppressor genes
TFTranscription factor
EMTEpithelial to mesenchymal transition
BELBiological expression language
GEOGene expression omnibus
LINCSLibrary of network-based integrated cellular signatures
NPANetwork perturbation amplitude
PPIProtein-protein interactions
KSKolmogorov-Smirnov test
GRmaxMaximum growth-rate inhibition
ECDFEmpirical cumulative distribution function

References

  1. Moffat, J.G.; Vincent, F.; Lee, J.A.; Eder, J.; Prunotto, M. Opportunities and challenges in phenotypic drug discovery: An industry perspective. Nat. Rev. Drug Discov. 2017, 16, 531–543. [Google Scholar] [CrossRef]
  2. Martin, F.; Thomson, T.M.; Sewer, A.; Drubin, D.A.; Mathis, C.; Weisensee, D.; Pratt, D.; Hoeng, J.; Peitsch, M.C. Assessment of network perturbation amplitudes by applying high-throughput data to causal biological networks. BMC Syst. Biol. 2012, 6, 1–18. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  3. Catlett, N.L.; Bargnesi, A.J.; Ungerer, S.; Seagaran, T.; Ladd, W.; Elliston, K.O.; Pratt, D. Reverse causal reasoning: Applying qualitative causal knowledge to the interpretation of high-throughput data. BMC Bioinform. 2013, 14, 340. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  4. Martin, F.; Sewer, A.; Talikka, M.; Xiang, Y.; Hoeng, J.; Peitsch, M.C. Quantification of biological network perturbations for mechanistic insight and diagnostics using two-layer causal models. BMC Bioinform. 2014, 15, 238. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  5. Barabasi, A.L.; Oltvai, Z.N. Network biology: Understanding the cell’s functional organization. Nat. Rev. Genet. 2004, 5, 101–113. [Google Scholar] [CrossRef]
  6. Shankavaram, U.T.; Reinhold, W.C.; Nishizuka, S.; Major, S.; Morita, D.; Chary, K.K.; Reimers, M.A.; Scherf, U.; Kahn, A.; Dolginow, D.; et al. Transcript and protein expression profiles of the NCI-60 cancer cell panel: An integromic microarray study. Mol. Cancer Ther. 2007, 6, 820–832. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  7. Ayala de la Peña, F.; Andrés, R.; Garcia-Sáenz, J.A.; Manso, L.; Margelí, M.; Dalmau, E.; Pernas, S.; Prat, A.; Servitja, S.; Ciruelos, E. SEOM clinical guidelines in early stage breast cancer (2018). Clin. Transl. Oncol. 2019, 21, 18–30. [Google Scholar] [CrossRef] [Green Version]
  8. Kogure, A.; Yoshioka, Y.; Ochiya, T. Extracellular vesicles in cancer metastasis: Potential as therapeutic targets and materials. Int. J. Mol. Sci. 2020, 21, 4463. [Google Scholar] [CrossRef] [PubMed]
  9. Khan, I.; Steeg, P.S. Metastasis suppressors: Functional pathways. Lab Investig. 2018, 98, 198–210. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  10. Keller, E.T. Metastasis suppressor genes: A role for raf kinase inhibitor protein (RKIP). Anticancer Drugs 2004, 15, 663–669. [Google Scholar] [CrossRef] [PubMed]
  11. Park, S.; Yeung, M.L.; Beach, S.; Shields, J.M.; Yeung, K.C. RKIP downregulates B-Raf kinase activity in melanoma cancer cells. Oncogene 2005, 24, 3535–3540. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  12. Zou, Q.; Wu, H.; Fu, F.; Yi, W.; Pei, L.; Zhou, M. RKIP suppresses the proliferation and metastasis of breast cancer cell lines through up-regulation of miR-185 targeting HMGA2. Arch. Biochem. Biophys. 2016, 610, 25–32. [Google Scholar] [CrossRef] [PubMed]
  13. Feng, L.; Zhang, C.; Liu, G.; Wang, F. RKIP negatively regulates the glucose induced angiogenesis and endothelial-mesenchymal transition in retinal endothelial cells. Exp. Eye Res. 2019, 189, 107851. [Google Scholar] [CrossRef] [PubMed]
  14. Chatterjee, D.; Bai, Y.; Wang, Z.; Beach, S.; Mott, S.; Roy, R.; Braastad, C.; Sun, Y.; Mukhopadhyay, A.; Aggarwal, B.B.; et al. RKIP Sensitizes Prostate and Breast Cancer Cells to Drug-induced Apoptosis. J. Biol. Chem. 2004, 279, 17515–17523. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  15. Szklarczyk, D.; Morris, J.H.; Cook, H.; Kuhn, M.; Wyder, S.; Simonovic, M.; Santos, A.; Doncheva, N.T.; Roth, A.; Bork, P.; et al. The STRING database in 2017: Quality-controlled protein-protein association networks, made broadly accessible. Nucleic Acids Res. 2017, 45, D362–D368. [Google Scholar] [CrossRef] [PubMed]
  16. Hoyt, C.T.; Domingo-Fernández, D.; Hofmann-Apitius, M. BEL Commons: An environment for exploration and analysis of networks encoded in Biological Expression Language. Database J. Biol. Databases Curation 2018, 2018, bay126. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  17. Yau, C.; Benz, C.C. Genes responsive to both oxidant stress and loss of estrogen receptor function identify a poor prognosis group of estrogen receptor positive primary breast cancers. Breast Cancer Res. BCR 2008, 10, R61. [Google Scholar] [CrossRef] [Green Version]
  18. Dahlman-Wright, K.; Qiao, Y.; Jonsson, P.; Gustafsson, J.A.; Williams, C.; Zhao, C. Interplay between AP-1 and estrogen receptor α in regulating gene expression and proliferation networks in breast cancer cells. Carcinogenesis 2012, 33, 1684–1691. [Google Scholar] [CrossRef] [Green Version]
  19. Bergamaschi, A.; Madak-Erdogan, Z.; Kim, Y.J.; Choi, Y.L.; Lu, H.; Katzenellenbogen, B.S. The forkhead transcription factor FOXM1 promotes endocrine resistance and invasiveness in estrogen receptor-positive breast cancer by expansion of stem-like cancer cells. Breast Cancer Res. BCR 2014, 16, 436. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  20. Theodorou, V.; Stark, R.; Menon, S.; Carroll, J.S. GATA3 acts upstream of FOXA1 in mediating ESR1 binding by shaping enhancer accessibility. Genome Res. 2013, 23, 12–22. [Google Scholar] [CrossRef] [Green Version]
  21. Elvidge, G.P.; Glenny, L.; Appelhoff, R.J.; Ratcliffe, P.J.; Ragoussis, J.; Gleadle, J.M. Concordant regulation of gene expression by hypoxia and 2-oxoglutarate-dependent dioxygenase inhibition: The role of HIF-1alpha, HIF-2alpha, and other pathways. J. Biol. Chem. 2006, 281, 15215–15226. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  22. Lai, C.F.; Flach, K.D.; Alexi, X.; Fox, S.P.; Ottaviani, S.; Thiruchelvam, P.T.R.; Kyle, F.J.; Thomas, R.S.; Launchbury, R.; Hua, H.; et al. Co-regulated gene expression by oestrogen receptor α and liver receptor homolog-1 is a feature of the oestrogen response in breast cancer cells. Nucleic Acids Res. 2013, 41, 10228–10240. [Google Scholar] [CrossRef] [PubMed]
  23. Salazar, M.D.; Ratnam, M.; Patki, M.; Kisovic, I.; Trumbly, R.; Iman, M.; Ratnam, M. During hormone depletion or tamoxifen treatment of breast cancer cells the estrogen receptor apoprotein supports cell cycling through the retinoic acid receptor α1 apoprotein. Breast Cancer Res. BCR 2011, 13, R18. [Google Scholar] [CrossRef] [Green Version]
  24. Buchwalter, G.; Hickey, M.M.; Cromer, A.; Selfors, L.M.; Gunawardane, R.N.; Frishman, J.; Jeselsohn, R.; Lim, E.; Chi, D.; Fu, X.; et al. PDEF promotes luminal differentiation and acts as a survival factor for ER-positive breast cancer cells. Cancer Cell 2013, 23, 753–767. [Google Scholar] [CrossRef] [Green Version]
  25. Tan, S.K.; Lin, Z.H.; Chang, C.W.; Varang, V.; Chng, K.R.; Pan, Y.F.; Yong, E.L.; Sung, W.K.; Sung, W.K.; Cheung, E. AP-2γ regulates oestrogen receptor-mediated long-range chromatin interaction and gene transcription. EMBO J. 2011, 30, 2569–2581. [Google Scholar] [CrossRef] [PubMed]
  26. Lasham, A.; Samuel, W.; Cao, H.; Patel, R.; Mehta, R.; Stern, J.L.; Reid, G.; Woolley, A.G.; Miller, L.D.; Black, M.A.; et al. YB-1, the E2F pathway, and regulation of tumor cell growth. J. Natl. Cancer Inst. 2012, 104, 133–146. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  27. Dunham, I.; Kundaje, A.; Aldred, S.F.; Collins, P.J.; Davis, C.A.; Doyle, F.; Epstein, C.B.; Frietze, S.; Harrow, J.; Kaul, R.; et al. An integrated encyclopedia of DNA elements in the human genome. Nature 2012, 489. [Google Scholar] [CrossRef]
  28. Marino, N.; Collins, J.W.; Shen, C.; Caplen, N.J.; Merchant, A.S.; Gökmen-Polar, Y.; Goswami, C.P.; Hoshino, T.; Qian, Y.; Sledge, G.W.; et al. Identification and validation of genes with expression patterns inverse to multiple metastasis suppressor genes in breast cancer cell lines. Clin. Exp. Metastasis 2014, 31, 771–786. [Google Scholar] [CrossRef] [PubMed]
  29. Sean, D.; Meltzer, P.S. GEOquery: A bridge between the Gene Expression Omnibus (GEO) and BioConductor. Bioinformatics 2007, 23, 1846–1847. [Google Scholar] [CrossRef] [Green Version]
  30. Ritchie, M.E.; Phipson, B.; Wu, D.; Hu, Y.; Law, C.W.; Shi, W.; Smyth, G.K. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015, 43, e47. [Google Scholar] [CrossRef] [PubMed]
  31. Feng, C.; Song, C.; Liu, Y.; Qian, F.; Gao, Y.; Ning, Z.; Wang, Q.; Jiang, Y.; Li, Y.; Li, M.; et al. KnockTF: A comprehensive human gene expression profile database with knockdown/knockout of transcription factors. Nucleic Acids Res. 2019, 48, D93–D100. [Google Scholar] [CrossRef]
  32. Hafner, M.; Niepel, M.; Chung, M.; Sorger, P.K. Growth rate inhibition metrics correct for confounders in measuring sensitivity to cancer drugs. Nat. Methods 2016, 13, 521–527. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  33. Koleti, A.; Terryn, R.; Stathias, V.; Chung, C.; Cooper, D.J.; Turner, J.P.; VidoviÄ, D.; Forlin, M.; Kelley, T.T.; D’Urso, A.; et al. Data Portal for the Library of Integrated Network-based Cellular Signatures (LINCS) program: Integrated access to diverse large-scale cellular perturbation response data. Nucleic Acids Res. 2018, 46, D558–D566. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  34. Kort, E. Slinky: Putting the Fun in LINCS L1000 Data Analysis. Available online: https://bioconductor.org/packages/release/bioc/html/slinky.html (accessed on 29 November 2021).
  35. Martin, F.; Gubian, S.; Talikka, M.; Hoeng, J.; Peitsch, M.C. NPA: An R package for computing network perturbation amplitudes using gene expression data and two-layer networks. BMC Bioinform. 2019, 20, 1–9. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  36. Baaten, B.J.; Li, C.R.; Bradley, L.M. Multifaceted regulation of T cells by CD44. Commun. Integr. Biol. 2010, 3, 508–512. [Google Scholar] [CrossRef] [Green Version]
  37. Wei, Q.; Zhang, F.; Richardson, M.M.; Roy, N.H.; Rodgers, W.; Liu, Y.; Zhao, W.; Fu, C.; Ding, Y.; Huang, C.; et al. CD82 restrains pathological angiogenesis by altering lipid raft clustering and CD44 trafficking in endothelial cells. Circulation 2014, 130, 1493–1504. [Google Scholar] [CrossRef] [Green Version]
  38. Tang, Y.; Cheng, Y.; Martinka, M.; Ong, C.J.; Li, G. Prognostic significance of KAI1/CD82 in human melanoma and its role in cell migration and invasion through the regulation of ING4. Carcinogenesis 2014, 35, 86–95. [Google Scholar] [CrossRef] [Green Version]
  39. Lee, M.S.; Byun, H.J.; Lee, J.; Jeoung, D.I.; Kim, Y.M.; Lee, H. Tetraspanin CD82 represses Sp1-mediated Snail expression and the resultant E-cadherin expression interrupts nuclear signaling of β-catenin by increasing its membrane localization. Cell. Signal. 2018, 52, 83–94. [Google Scholar] [CrossRef] [PubMed]
  40. Lu, M.; Marsters, S.; Ye, X.; Luis, E.; Gonzalez, L.; Ashkenazi, A. E-cadherin couples death receptors to the cytoskeleton to regulate apoptosis. Mol. Cell 2014, 54, 987–998. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  41. Kim, J.C.; Ha, Y.J.; Tak, K.H.; Roh, S.A.; Kwon, Y.H.; Kim, C.W.; Yoon, Y.S.; Lee, J.L.; Park, Y.; Kim, S.K.; et al. Opposite functions of GSN and OAS2 on colorectal cancer metastasis, mediating perineural and lymphovascular invasion, respectively. PLoS ONE 2018, 13, e0202856. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  42. Raingeaud, J.; Whitmarsh, A.J.; Barrett, T.; Dérijard, B.; Davis, R.J. MKK3- and MKK6-regulated gene expression is mediated by the p38 mitogen-activated protein kinase signal transduction pathway. Mol. Cell. Biol. 1996, 16, 1247–1255. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  43. Haeusgen, W.; Tueffers, L.; Herdegen, T.; Waetzig, V. Map2k4δ-identification and functional characterization of a novel Map2k4 splice variant. Biochim. Biophys. Acta 2014, 1843, 875–884. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  44. El Mchichi, B.; Hadji, A.; Vazquez, A.; Leca, G. p38 MAPK and MSK1 mediate caspase-8 activation in manganese-induced mitochondria-dependent cell death. Cell Death Differ. 2007, 14, 1826–1836. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  45. Shintani, Y.; Fukumoto, Y.; Chaika, N.; Svoboda, R.; Wheelock, M.J.; Johnson, K.R. Collagen I-mediated up-regulation of N-cadherin requires cooperative signals from integrins and discoidin domain receptor 1. J. Cell Biol. 2008, 180, 1277–1289. [Google Scholar] [CrossRef] [Green Version]
  46. Ventura, J.J.; Tenbaum, S.; Perdiguero, E.; Huth, M.; Guerra, C.; Barbacid, M.; Pasparakis, M.; Nebreda, A.R. p38alpha MAP kinase is essential in lung stem and progenitor cell proliferation and differentiation. Nat. Genet. 2007, 39, 750–758. [Google Scholar] [CrossRef] [PubMed]
  47. Alvarado-Kristensson, M.; Melander, F.; Leandersson, K.; Rönnstrand, L.; Wernstedt, C.; Andersson, T. p38-MAPK signals survival by phosphorylation of caspase-8 and caspase-3 in human neutrophils. J. Exp. Med. 2004, 199, 449–458. [Google Scholar] [CrossRef] [PubMed]
  48. Janknecht, R.; Hunter, T. Convergence of MAP kinase pathways on the ternary complex factor Sap-1a. EMBO J. 1997, 16, 1620–1627. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  49. Hutchison, M.R. Mice with a conditional deletion of the neurotrophin receptor TrkB are dwarfed, and are similar to mice with a MAPK14 deletion. PLoS ONE 2013, 8, e66206. [Google Scholar] [CrossRef]
  50. McCorkle, J.R.; Leonard, M.K.; Kraner, S.D.; Blalock, E.M.; Ma, D.; Zimmer, S.G.; Kaetzel, D.M. The metastasis suppressor NME1 regulates expression of genes linked to metastasis and patient outcome in melanoma and breast carcinoma. Cancer Genom. Proteom. 2014, 11, 175–194. [Google Scholar]
  51. Berger, J.C.; Vander Griend, D.J.; Robinson, V.L.; Hickson, J.A.; Rinker-Schaeffer, C.W. Metastasis suppressor genes: From gene identification to protein function and regulation. Cancer Biol. Ther. 2005, 4, 805–812. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  52. Lai, R.; Gu, M.; Jiang, W.; Lin, W.; Xu, P.; Liu, Z.; Huang, H.; An, H.; Wang, X. Raf Kinase Inhibitor Protein Preferentially Promotes TLR3-Triggered Signaling and Inflammation. J. Immunol. 2017, 198, 4086–4095. [Google Scholar] [CrossRef] [Green Version]
  53. Wu, Y.; Chen, K.; Xing, G.; Li, L.; Ma, B.; Hu, Z.; Duan, L.; Liu, X. Phospholipid remodeling is critical for stem cell pluripotency by facilitating mesenchymal-to-epithelial transition. Sci. Adv. 2019, 5, eaax7525. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  54. Saito, K.; Takahashi, K.; Huang, B.; Asahara, M.; Kiso, H.; Togo, Y.; Tsukamoto, H.; Mishima, S.; Nagata, M.; Iida, M.; et al. Loss of Stemness, EMT, and Supernumerary Tooth Formation in Cebpb-/-Runx2+/− Murine Incisors. Sci. Rep. 2018, 8, 5169. [Google Scholar] [CrossRef] [PubMed]
  55. Mishra, A.K.; Parish, C.R.; Wong, M.L.; Licinio, J.; Blackburn, A.C. Leptin signals via TGFB1 to promote metastatic potential and stemness in breast cancer. PLoS ONE 2017, 12, e0178454. [Google Scholar] [CrossRef] [Green Version]
  56. Wang, H.; Leinwand, L.A.; Anseth, K.S. Roles of transforming growth factor-β1 and OB-cadherin in porcine cardiac valve myofibroblast differentiation. Faseb J. Off. Publ. Fed. Am. Soc. Exp. Biol. 2014, 28, 4551–4562. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  57. Liu, X.; Su, L.; Liu, X. Loss of CDH1 up-regulates epidermal growth factor receptor via phosphorylation of YBX1 in non-small cell lung cancer cells. FEBS Lett. 2013, 587, 3995–4000. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  58. Chen, Z.; Li, L.; Xu, S.; Liu, Z.; Zhou, C.; Li, Z.; Liu, Y.; Wu, W.; Huang, Y.; Kuang, M.; et al. A Cdh1-FoxM1-Apc axis controls muscle development and regeneration. Cell Death Dis. 2020, 11, 180. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  59. Fujita, N.; Kajita, M.; Taysavang, P.; Wade, P.A. Hormonal regulation of metastasis-associated protein 3 transcription in breast cancer cells. Mol. Endocrinol. 2004, 18, 2937–2949. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  60. Cardamone, M.D.; Bardella, C.; Gutierrez, A.; Di Croce, L.; Rosenfeld, M.G.; Di Renzo, M.F.; De Bortoli, M. ERalpha as ligand-independent activator of CDH-1 regulates determination and maintenance of epithelial morphology in breast cancer cells. Proc. Natl. Acad. Sci. USA 2009, 106, 7420–7425. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  61. Ma, H.; Gollahon, L.S. ERα Mediates Estrogen-Induced Expression of the Breast Cancer Metastasis Suppressor Gene BRMS1. Int. J. Mol. Sci. 2016, 17, 158. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  62. Fujita, N.; Jaye, D.L.; Kajita, M.; Geigerman, C.; Moreno, C.S.; Wade, P.A. MTA3, a Mi-2/NuRD complex subunit, regulates an invasive growth pathway in breast cancer. Cell 2003, 113, 207–219. [Google Scholar] [CrossRef] [Green Version]
  63. Kammerer, M.; Gutzwiller, S.; Stauffer, D.; Delhon, I.; Seltenmeyer, Y.; Fournier, B. Estrogen Receptor α (ERα) and Estrogen Related Receptor α (ERRα) are both transcriptional regulators of the Runx2-I isoform. Mol. Cell. Endocrinol. 2013, 369, 150–160. [Google Scholar] [CrossRef] [PubMed]
  64. Hovhannisyan, H.; Zhang, Y.; Hassan, M.Q.; Wu, H.; Glackin, C.; Lian, J.B.; Stein, J.L.; Montecino, M.; Stein, G.S.; van Wijnen, A.J. Genomic occupancy of HLH, AP1 and Runx2 motifs within a nuclease sensitive site of the Runx2 gene. J. Cell. Physiol. 2013, 228, 313–321. [Google Scholar] [CrossRef] [Green Version]
  65. Dong, C.; Ye, D.X.; Zhang, W.B.; Pan, H.Y.; Zhang, Z.Y.; Zhang, L. Overexpression of c-fos promotes cell invasion and migration via CD44 pathway in oral squamous cell carcinoma. J. Oral Pathol. Med. Off. Publ. Int. Assoc. Oral Pathol. Am. Acad. Oral Pathol. 2015, 44, 353–360. [Google Scholar] [CrossRef]
  66. Chen, H.J.; Huang, R.L.; Liew, P.L.; Su, P.H.; Chen, L.Y.; Weng, Y.C.; Chang, C.C.; Wang, Y.C.; Chan, M.W.Y.; Lai, H.C. GATA3 as a master regulator and therapeutic target in ovarian high-grade serous carcinoma stem cells. Int. J. Cancer 2018, 143, 3106–3119. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  67. Krishnamachary, B.; Penet, M.F.; Nimmagadda, S.; Mironchik, Y.; Raman, V.; Solaiyappan, M.; Semenza, G.L.; Pomper, M.G.; Bhujwalla, Z.M. Hypoxia regulates CD44 and its variant isoforms through HIF-1α in triple negative breast cancer. PLoS ONE 2012, 7, e44078. [Google Scholar] [CrossRef]
  68. Xue, Z.; Vis, D.J.; Bruna, A.; Sustic, T.; van Wageningen, S.; Batra, A.S.; Rueda, O.M.; Bosdriesz, E.; Caldas, C.; Wessels, L.F.A.; et al. MAP3K1 and MAP2K4 mutations are associated with sensitivity to MEK inhibitors in multiple cancer models. Cell Res. 2018, 28, 719–729. [Google Scholar] [CrossRef]
  69. Wan, Y.Y. GATA3: A master of many trades in immune regulation. Trends Immunol. 2014, 35, 233–242. [Google Scholar] [CrossRef] [PubMed]
  70. Liu, Y.; Mayo, M.W.; Nagji, A.S.; Smith, P.W.; Ramsey, C.S.; Li, D.; Jones, D.R. Phosphorylation of RelA/p65 promotes DNMT-1 recruitment to chromatin and represses transcription of the tumor metastasis suppressor gene BRMS1. Oncogene 2012, 31, 1143–1154. [Google Scholar] [CrossRef]
  71. Lee, S.H.; Che, X.; Jeong, J.H.; Choi, J.Y.; Lee, Y.J.; Lee, Y.H.; Bae, S.C.; Lee, Y.M. Runx2 protein stabilizes hypoxia-inducible factor-1α through competition with von Hippel-Lindau protein (pVHL) and stimulates angiogenesis in growth plate hypertrophic chondrocytes. J. Biol. Chem. 2012, 287, 14760–14771. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  72. Shi, N.; Zhang, J.; Chen, S.Y. Runx2, a novel regulator for goblet cell differentiation and asthma development. Faseb J. Off. Publ. Fed. Am. Soc. Exp. Biol. 2017, 31, 412–420. [Google Scholar] [CrossRef] [Green Version]
  73. Han, H.J.; Russo, J.; Kohwi, Y.; Kohwi-Shigematsu, T. SATB1 reprogrammes gene expression to promote breast tumour growth and metastasis. Nature 2008, 452, 187–193. [Google Scholar] [CrossRef] [PubMed]
  74. Very, N.; Lefebvre, T.; El Yazidi-Belkoura, I. Drug resistance related to aberrant glycosylation in colorectal cancer. Oncotarget 2018, 9, 1380–1402. [Google Scholar] [CrossRef] [Green Version]
  75. Beach, S.; Tang, H.; Park, S.; Dhillon, A.S.; Keller, E.T.; Kolch, W.; Yeung, K.C. Snail is a repressor of RKIP transcription in metastatic prostate cancer cells. Oncogene 2008, 27, 2243–2248. [Google Scholar] [CrossRef] [Green Version]
  76. Vega, S.; Morales, A.V.; Ocaña, O.H.; Valdés, F.; Fabregat, I.; Nieto, M.A. Snail blocks the cell cycle and confers resistance to cell death. Genes Dev. 2004, 18, 1131–1143. [Google Scholar] [CrossRef] [Green Version]
  77. Spanheimer, P.M.; Askeland, R.W.; Kulak, M.V.; Wu, T.; Weigel, R.J. High TFAP2C/low CD44 expression is associated with an increased rate of pathologic complete response following neoadjuvant chemotherapy in breast cancer. J. Surg. Res. 2013, 184, 519–525. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  78. Hsieh, A.C.; Liu, Y.; Edlind, M.P.; Ingolia, N.T.; Janes, M.R.; Sher, A.; Shi, E.Y.; Stumpf, C.R.; Christensen, C.; Bonham, M.J.; et al. The translational landscape of mTOR signalling steers cancer initiation and metastasis. Nature 2012, 485, 55–61. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  79. Yang, J.; AlTahan, A.; Jones, D.T.; Buffa, F.M.; Bridges, E.; Interiano, R.B.; Qu, C.; Vogt, N.; Li, J.L.; Baban, D.; et al. Estrogen receptor-α directly regulates the hypoxia-inducible factor 1 pathway associated with antiestrogen response in breast cancer. Proc. Natl. Acad. Sci. USA 2015, 112, 15172–15177. [Google Scholar] [CrossRef] [Green Version]
  80. Carr, J.R.; Kiefer, M.M.; Park, H.J.; Li, J.; Wang, Z.; Fontanarosa, J.; DeWaal, D.; Kopanja, D.; Benevolenskaya, E.V.; Guzman, G.; et al. FoxM1 regulates mammary luminal cell fate. Cell Rep. 2012, 1, 715–729. [Google Scholar] [CrossRef] [Green Version]
  81. Du, L.; Ning, Z.; Zhang, H.; Liu, F. Corepressor metastasis-associated protein 3 modulates epithelial-to-mesenchymal transition and metastasis. Chin. J. Cancer 2017, 36, 28. [Google Scholar] [CrossRef] [Green Version]
  82. Thiruchelvam, P.T.R.; Lai, C.F.; Hua, H.; Thomas, R.S.; Hurtado, A.; Hudson, W.; Bayly, A.R.; Kyle, F.J.; Periyasamy, M.; Photiou, A.; et al. The liver receptor homolog-1 regulates estrogen receptor expression in breast cancer cells. Breast Cancer Res. Treat. 2011, 127, 385–396. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  83. Pfahl, M. Nuclear receptor/AP-1 interaction. Endocr. Rev. 1993, 14, 651–658. [Google Scholar] [CrossRef] [PubMed]
  84. Julien, S.; Puig, I.; Caretti, E.; Bonaventure, J.; Nelles, L.; van Roy, F.; Dargemont, C.; de Herreros, A.G.; Bellacosa, A.; Larue, L. Activation of NF-kappaB by Akt upregulates Snail expression and induces epithelium mesenchyme transition. Oncogene 2007, 26, 7445–7456. [Google Scholar] [CrossRef] [PubMed]
  85. Woodfield, G.W.; Chen, Y.; Bair, T.B.; Domann, F.E.; Weigel, R.J. Identification of primary gene targets of TFAP2C in hormone responsive breast carcinoma cells. Genes Chromosom. Cancer 2010, 49, 948–962. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  86. Mizamtsidi, M.; Nastos, C.; Mastorakos, G.; Dina, R.; Vassiliou, I.; Gazouli, M.; Palazzo, F. Diagnosis, management, histology and genetics of sporadic primary hyperparathyroidism: Old knowledge with new tricks. Endocr. Connect. 2018, 7, R56–R68. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  87. Ranganathan, P.; Pramesh, C.; Aggarwal, R. Common pitfalls in statistical analysis: Measures of agreement. Perspect. Clin. Res. 2017, 8. [Google Scholar] [CrossRef]
  88. Ahmed, M.; Kim, D.R. pcr: An R package for quality assessment, analysis and testing of qPCR data. PeerJ 2018, 6, e4473. [Google Scholar] [CrossRef]
  89. R Core Team. R: A Language and Environment for Statistical Computing; R Core Team: Vienna, Austria, 2021. [Google Scholar]
  90. Huber, W.; Carey, V.J.; Gentleman, R.; Anders, S.; Carlson, M.; Carvalho, B.S.; Bravo, H.C.; Davis, S.; Gatto, L.; Girke, T.; et al. Orchestrating high-throughput genomic analysis with Bioconductor. Nat. Methods 2015, 12, 115–121. [Google Scholar] [CrossRef]
  91. Ruiz Castro, P.A.; Yepiskoposyan, H.; Gubian, S.; Calvino-Martin, F.; Kogel, U.; Renggli, K.; Peitsch, M.C.; Hoeng, J.; Talikka, M. Systems biology approach highlights mechanistic differences between Crohn’s disease and ulcerative colitis. Sci. Rep. 2021, 11, 11519. [Google Scholar] [CrossRef]
  92. Lo Sasso, G.; Phillips, B.W.; Sewer, A.; Battey, J.N.; Kondylis, A.; Talikka, M.; Titz, B.; Guedj, E.; Peric, D.; Bornand, D.; et al. The reduction of DSS-induced colitis severity in mice exposed to cigarette smoke is linked to immune modulation and microbial shifts. Sci. Rep. 2020, 10, 3829. [Google Scholar] [CrossRef] [Green Version]
  93. Martin, T.A.; Goyal, A.; Watkins, G.; Jiang, W.G. Expression of the transcription factors snail, slug, and twist and their clinical significance in human breast cancer. Ann. Surg. Oncol. 2005, 12, 488–496. [Google Scholar] [CrossRef] [PubMed]
  94. Barberà, M.J.; Puig, I.; Domínguez, D.; Julien-Grille, S.; Guaita-Esteruelas, S.; Peiró, S.; Baulida, J.; Francí, C.; Dedhar, S.; Larue, L.; et al. Regulation of Snail transcription during epithelial to mesenchymal transition of tumor cells. Oncogene 2004, 23, 7345–7354. [Google Scholar] [CrossRef] [Green Version]
  95. Wu, K.; Bonavida, B. The activated NF-kappaB-Snail-RKIP circuitry in cancer regulates both the metastatic cascade and resistance to apoptosis by cytotoxic drugs. Crit. Rev. Immunol. 2009, 29, 241–254. [Google Scholar] [CrossRef] [PubMed]
  96. You, D.J.; Park, C.R.; Lee, H.B.; Moon, M.J.; Kang, J.H.; Lee, C.; Oh, S.H.; Ahn, C.; Seong, J.Y.; Hwang, J.I. A splicing variant of NME1 negatively regulates NF-κB signaling and inhibits cancer metastasis by interacting with IKKβ. J. Biol. Chem. 2014, 289, 17709–17720. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Cancers 13 06098 g001 550
Figure 1. A workflow for building, refining, and evaluating a breast cancer-specific network model of metastasis. (A) We queried the String database for protein–protein interactions for metastasis suppressor genes (MSG) and transcription factors (TF). Next, we surveyed the relevant literature to vet and supplement the interactions and express them as possible causal links in the biological expression language (BEL). Finally, we filtered and evaluated the network model using knockdown and drug perturbation datasets. (B) Conflicts between the data-driven and curated interactions were resolved by prioritizing the context, the strength of the evidence, and the effect sizes. We scored the metastasis model on drug perturbations using the network perturbation amplitude (NPA). We divided the network into smaller subgraphs and measured the agreement between the expected and the observed direction activity. We considered the consistency in (C) subnetworks of nodes connected to an upstream by one edge (concordance) and (D) paths connecting a particular node to its upstream by a sequence of edges (coherence).
Figure 1. A workflow for building, refining, and evaluating a breast cancer-specific network model of metastasis. (A) We queried the String database for protein–protein interactions for metastasis suppressor genes (MSG) and transcription factors (TF). Next, we surveyed the relevant literature to vet and supplement the interactions and express them as possible causal links in the biological expression language (BEL). Finally, we filtered and evaluated the network model using knockdown and drug perturbation datasets. (B) Conflicts between the data-driven and curated interactions were resolved by prioritizing the context, the strength of the evidence, and the effect sizes. We scored the metastasis model on drug perturbations using the network perturbation amplitude (NPA). We divided the network into smaller subgraphs and measured the agreement between the expected and the observed direction activity. We considered the consistency in (C) subnetworks of nodes connected to an upstream by one edge (concordance) and (D) paths connecting a particular node to its upstream by a sequence of edges (coherence).
Cancers 13 06098 g001
Cancers 13 06098 g002 550
Figure 2. Fold-change between control and knockdown of metastasis suppressors and transcription factors. We obtained expression profiles of seventeen metastasis suppressor genes (MSG) and eleven transcription factor (TF) knockdowns in the MCF7 cell line. We calculated the differential expression between the knockdown and control samples. Bottom, bars represent the fold-change (log 2 ) of the target gene (Self) between the control and knockdown. Top, box plots represent the distribution of the fold-change of (Other) MSGs and TFs in the corresponding knockdown condition.
Figure 2. Fold-change between control and knockdown of metastasis suppressors and transcription factors. We obtained expression profiles of seventeen metastasis suppressor genes (MSG) and eleven transcription factor (TF) knockdowns in the MCF7 cell line. We calculated the differential expression between the knockdown and control samples. Bottom, bars represent the fold-change (log 2 ) of the target gene (Self) between the control and knockdown. Top, box plots represent the distribution of the fold-change of (Other) MSGs and TFs in the corresponding knockdown condition.
Cancers 13 06098 g002
Cancers 13 06098 g003 550
Figure 3. Network of interactions among metastasis suppressors and their regulators. (A) A network graph of seventeen metastasis suppressors and eleven transcription factor nodes. We curated the edges (red, repression or blue, activation) from protein–protein interactions, literature survey, and context-specific gene expression data (Functional layer). (B) Top left, the fractions of nodes with a given number of edges (degree). Top right, the density function of the sum of the shortest paths between every node and all others (closeness). We built a transcript layer of significant changes in gene expression (absolute log 2 fold-change > 0.5 and p-value < 0.01) as a result of knocking down every node in the functional layer. Bottom left, the density function of the numbers of nodes in the transcript layer connected to the nodes in the network. Bottom right, a histogram of the numbers up-and down-regulated nodes in the transcript layer.
Figure 3. Network of interactions among metastasis suppressors and their regulators. (A) A network graph of seventeen metastasis suppressors and eleven transcription factor nodes. We curated the edges (red, repression or blue, activation) from protein–protein interactions, literature survey, and context-specific gene expression data (Functional layer). (B) Top left, the fractions of nodes with a given number of edges (degree). Top right, the density function of the sum of the shortest paths between every node and all others (closeness). We built a transcript layer of significant changes in gene expression (absolute log 2 fold-change > 0.5 and p-value < 0.01) as a result of knocking down every node in the functional layer. Bottom left, the density function of the numbers of nodes in the transcript layer connected to the nodes in the network. Bottom right, a histogram of the numbers up-and down-regulated nodes in the transcript layer.
Cancers 13 06098 g003
Cancers 13 06098 g004 550
Figure 4. Growth inhibition of MFC7 and metastasis network perturbation under drug treatments. We obtained growth inhibition and gene expression of breast cancer cell lines (n = 71) treated with different drugs (n = 35). (A) Maximum growth rate inhibition (GRmax) of MCF7 from different does and replicates after accounting for the baseline growth rate. Values of 1 indicate no inhibition, while −1 indicates the maximum inhibition of growth. (B) A histogram of the Pearson’s correlation coefficients between MCF7 GRmax and other cell lines. (C) We computed the network perturbation amplitudes (NPA) for the metastasis network model and every node in the network under drug (n = 35) treatment in different cell lines (n = 67). Positive values indicate the drug treatment activates the node in MCF7 and negative values indicate repression of the node. (D) A histogram of the Pearson’s correlation coefficients between MCF7 NPA and other cell lines.
Figure 4. Growth inhibition of MFC7 and metastasis network perturbation under drug treatments. We obtained growth inhibition and gene expression of breast cancer cell lines (n = 71) treated with different drugs (n = 35). (A) Maximum growth rate inhibition (GRmax) of MCF7 from different does and replicates after accounting for the baseline growth rate. Values of 1 indicate no inhibition, while −1 indicates the maximum inhibition of growth. (B) A histogram of the Pearson’s correlation coefficients between MCF7 GRmax and other cell lines. (C) We computed the network perturbation amplitudes (NPA) for the metastasis network model and every node in the network under drug (n = 35) treatment in different cell lines (n = 67). Positive values indicate the drug treatment activates the node in MCF7 and negative values indicate repression of the node. (D) A histogram of the Pearson’s correlation coefficients between MCF7 NPA and other cell lines.
Cancers 13 06098 g004
Cancers 13 06098 g005 550
Figure 5. Concordance of expectations and observations in the subnetworks of the metastasis model. We calculated and binarized (1 for activation or −1 for repression) the perturbation coefficients of every node in the network. We then evaluated the agreement between the expected and observed direction of change in the subnetworks of nodes connected to an upstream by one edge (concordance). First, we multiplied the drug’s effect on the upstream node by the sign of the edges to form expectations. Next, we compared the expectations with the actual perturbation coefficients of the corresponding nodes. Negative Cohen’s κ indicates worse and positive better agreement than expected by chance. Finally, we compared the probability distribution of the concordance to randomly generated values using Kolmogorov-Smirnov (KS) test. D + is the maximum distance between the cumulative distribution functions (ECDF). (A) A histogram of the average concordances for every drug. (B) A histogram of the average concordances for every node. (C) A similar workflow was applied to other cancer cell lines (n = 67), and the concordance values averaged by cell lines are shown as a histogram.
Figure 5. Concordance of expectations and observations in the subnetworks of the metastasis model. We calculated and binarized (1 for activation or −1 for repression) the perturbation coefficients of every node in the network. We then evaluated the agreement between the expected and observed direction of change in the subnetworks of nodes connected to an upstream by one edge (concordance). First, we multiplied the drug’s effect on the upstream node by the sign of the edges to form expectations. Next, we compared the expectations with the actual perturbation coefficients of the corresponding nodes. Negative Cohen’s κ indicates worse and positive better agreement than expected by chance. Finally, we compared the probability distribution of the concordance to randomly generated values using Kolmogorov-Smirnov (KS) test. D + is the maximum distance between the cumulative distribution functions (ECDF). (A) A histogram of the average concordances for every drug. (B) A histogram of the average concordances for every node. (C) A similar workflow was applied to other cancer cell lines (n = 67), and the concordance values averaged by cell lines are shown as a histogram.
Cancers 13 06098 g005
Cancers 13 06098 g006 550
Figure 6. Coherence of expectations and observations in the paths to PEBP1. We calculated and binarized (1 for activation or −1 for repression) the perturbation coefficients of every node in the network. We then evaluated the agreement between the expected and observed direction of change in the paths connecting a PEBP1 to its upstream by a sequence of edges (coherence). First, we multiplied the drug’s effect on the upstream nodes by the sign of the edge connecting it to the next node to form expectations. Next, we compared the expectations with the actual perturbation coefficients of the corresponding nodes. Negative Cohen’s κ indicates worse and positive better agreement than expected by chance. Finally, we compared the probability distribution of the coherence to randomly generated values using Kolmogorov–Smirnov (KS) test. D + is the maximum distance between the cumulative distribution functions (ECDF). (A) A histogram of the average coherence for every drug. (B) A histogram of the average coherence for every path. (C) A similar workflow was applied to other cancer cell lines (n = 67), and the coherence values averaged by cell line are shown as a histogram.
Figure 6. Coherence of expectations and observations in the paths to PEBP1. We calculated and binarized (1 for activation or −1 for repression) the perturbation coefficients of every node in the network. We then evaluated the agreement between the expected and observed direction of change in the paths connecting a PEBP1 to its upstream by a sequence of edges (coherence). First, we multiplied the drug’s effect on the upstream nodes by the sign of the edge connecting it to the next node to form expectations. Next, we compared the expectations with the actual perturbation coefficients of the corresponding nodes. Negative Cohen’s κ indicates worse and positive better agreement than expected by chance. Finally, we compared the probability distribution of the coherence to randomly generated values using Kolmogorov–Smirnov (KS) test. D + is the maximum distance between the cumulative distribution functions (ECDF). (A) A histogram of the average coherence for every drug. (B) A histogram of the average coherence for every path. (C) A similar workflow was applied to other cancer cell lines (n = 67), and the coherence values averaged by cell line are shown as a histogram.
Cancers 13 06098 g006
Cancers 13 06098 g007 550
Figure 7. A functional model of PEBP1 interactions with other metastasis suppressors and regulators. We calculated the coherence of the paths leading to PEBP1 in the metastasis networks the expected to the observed perturbations with various drug treatments. (A) Top five paths based on the average coherence across drug treatments. (B) A network graph of the top coherent PEBP1 paths. Nodes represent biological entities and connect to each other by activation (arrow) or repression (blunt) edges. The thickness of the edges indicates the coherence of the path averaged by drugs.
Figure 7. A functional model of PEBP1 interactions with other metastasis suppressors and regulators. We calculated the coherence of the paths leading to PEBP1 in the metastasis networks the expected to the observed perturbations with various drug treatments. (A) Top five paths based on the average coherence across drug treatments. (B) A network graph of the top coherent PEBP1 paths. Nodes represent biological entities and connect to each other by activation (arrow) or repression (blunt) edges. The thickness of the edges indicates the coherence of the path averaged by drugs.
Cancers 13 06098 g007
Cancers 13 06098 g008 550
Figure 8. Activation and repression of PEBP1/RKIP and its regulatory pathways. MCF7 cells (n > 5) were treated with three activators (Epirubicin, Methotrexate, or Vorinostat), three repressors (Cisplatin, Imatinib, or Sorafenib), or the control DMSO for 24 h at 1.5 μ M dose. Total RNA was extracted and amplified using PCR. The amount of five mRNA (RELA, SNAI1, PEBP1, NME1, and ESR1) was quantified relative to the control GAPDH using RT-qPCR. (A) The relative RKIP mRNA with all six drug treatments compared to the control. The relative amount of five mRNA in three repressors treatments (B) Sorafenib, (C) Cisplatin or (D) Imatinib, and one activator (E) Epirubicin. (F) A model of the pathways regulating RKIP and the potential activators and repressors. Drugs are positioned near the node they regulate with plus or minus signs to indicate activation or repression. * indicates p < 0.05.
Figure 8. Activation and repression of PEBP1/RKIP and its regulatory pathways. MCF7 cells (n > 5) were treated with three activators (Epirubicin, Methotrexate, or Vorinostat), three repressors (Cisplatin, Imatinib, or Sorafenib), or the control DMSO for 24 h at 1.5 μ M dose. Total RNA was extracted and amplified using PCR. The amount of five mRNA (RELA, SNAI1, PEBP1, NME1, and ESR1) was quantified relative to the control GAPDH using RT-qPCR. (A) The relative RKIP mRNA with all six drug treatments compared to the control. The relative amount of five mRNA in three repressors treatments (B) Sorafenib, (C) Cisplatin or (D) Imatinib, and one activator (E) Epirubicin. (F) A model of the pathways regulating RKIP and the potential activators and repressors. Drugs are positioned near the node they regulate with plus or minus signs to indicate activation or repression. * indicates p < 0.05.
Cancers 13 06098 g008
Table
Table 1. Metastasis suppressor genes.
Table 1. Metastasis suppressor genes.
CategoryGenes
Cell-cell adhesionEpican (CD44), Tetraspanin 27 (CD82), Cadherin 11 (CDH11), Cadherin 2 (CDH2), Cadherin 1 (CDH1) and Gelsolin (GSN)
ScaffoldingGravin/a-kinase anchor protein 12 (AKAP12)
MAPKDual specificity mitogen-activated protein kinase kinase 6 (MAP2K6), 4 (MAP2K4), 7 (MAP2K7) and Mitogen-activated protein kinase 14 (MAPK14)
TranscriptionNME/NM23 Nucleoside Diphosphate Kinase 1 (NME1) and breast cancer metastasis-suppressor (BRMS1)
GTP-bindingRho GDP Dissociation Inhibitor Beta (ARGHDIB) and Developmentally-regulated GTP-biding protein 1 (DRG1)
OtherRibonucleotide Reductase Catalytic Subunit M1 (RRM1) and Phosphatidylethanolamine-binding protein 1 (PEBP1)
Table
Table 2. Transcription factors targeting metastasis suppressor genes in MCF7.
Table 2. Transcription factors targeting metastasis suppressor genes in MCF7.
TFNameDataset IDRef.
ESR1Estrogen receptor 1GSE10061[17]
FOSFos Proto-Oncogene AP-1 Transcription Factor SubunitGSE36586[18]
FOXM1Forkhead Box M1GSE55204[19]
GATA3GATA Binding Protein 3GSE39623[20]
HIF1AHypoxia Inducible Factor 1 Subunit AlphaGSE3188[21]
NR5A2Nuclear Receptor Subfamily 5 Group A Member 2GSE47803[22]
RARARetinoic Acid Receptor AlphaGSE26298[23]
SPDEFSAM Pointed Domain Containing ETS Transcription FactorGSE40985[24]
TFAP2CTranscription Factor AP-2 GammaGSE26740[25]
YBX1Y-Box Binding Protein 1GSE28433[26]
ZFXZinc Finger Protein X-LinkedENCSR005AHI[27]
Table
Table 4. RT-qPCR primers.
Table 4. RT-qPCR primers.
GeneForward PrimerReverse Primer
ESR15′-TGGAGTCTGGTCCTGTGAGG-3′5′-GGTCTTTTCGTATCCCACCTTTC-3′
SNAI15′-CCAGTGCCTCGACCACTATG-3′5′-CTGCTGGAAGGTAAACTCTGG-3′
RELA5′-CCTATAGAAGAGCAGCGTGGG-3′5′-AGATCTTGAGCTCGGCAGTG-3′
NME15′-ACTAAGTCAGCCTGGTGTGC-3′5′-CGCCTTGAAAGACGATCCCT-3′
PEBP15′-GTCACACTTTAGCGGCCTGT-3′5′-CTCTCCGATTATGTGGGCTC-3′
GAPDH5′-TGCACCACCAACTGCTTAGC-3′5′-GGCATGGACTGTGGTCATGAG-3′
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Ahmed, M.; Lai, T.H.; Kim, W.; Kim, D.R. A Functional Network Model of the Metastasis Suppressor PEBP1/RKIP and Its Regulators in Breast Cancer Cells. Cancers 2021, 13, 6098. https://doi.org/10.3390/cancers13236098

AMA Style

Ahmed M, Lai TH, Kim W, Kim DR. A Functional Network Model of the Metastasis Suppressor PEBP1/RKIP and Its Regulators in Breast Cancer Cells. Cancers. 2021; 13(23):6098. https://doi.org/10.3390/cancers13236098

Chicago/Turabian Style

Ahmed, Mahmoud, Trang Huyen Lai, Wanil Kim, and Deok Ryong Kim. 2021. "A Functional Network Model of the Metastasis Suppressor PEBP1/RKIP and Its Regulators in Breast Cancer Cells" Cancers 13, no. 23: 6098. https://doi.org/10.3390/cancers13236098

APA Style

Ahmed, M., Lai, T. H., Kim, W., & Kim, D. R. (2021). A Functional Network Model of the Metastasis Suppressor PEBP1/RKIP and Its Regulators in Breast Cancer Cells. Cancers, 13(23), 6098. https://doi.org/10.3390/cancers13236098

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop