{'fraction_non_alphanumeric': 0.02146715951094683, 'fraction_numerical': 0.008707705430764856, 'mean_word_length': 5.269385026737968, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Potential Breast Anticancer Drug Targets Revealed by Differential Gene Regulatory Network Analysis and Molecular Docking: Neoadjuvant Docetaxel Drug as a Case Study

Understanding gene-gene interaction and its causal relationship to protein-protein interaction is a viable route for understanding drug action at the genetic level, which is largely hindered by inability to robustly map gene regulatory networks. Here, we use biological prior knowledge of family-to-family gene interactions available in the KEGG database to reveal individual gene-to-gene interaction networks that underlie the gene expression profiles of 2 cell line data sets, sensitive and resistive to neoadjuvant docetaxel breast anticancer drug. Comparison of the topology of the 2 networks revealed that the resistant network is highly connected with 2 large domains of connectivity: one in which the RAF1 and MAP2K2 genes form hubs of connectivity and another in which the RAS gene is highly connected. On the contrary, the sensitive network is highly disrupted with a lower degree of connectivity. We investigated the interactions of the neoadjuvant docetaxel drug with the protein chains encoded by gene-gene interactions that underlie the disruption of the sensitive network topology using protein-protein and drug-protein docking techniques. We found that the sensitive network is likely to be disrupted by interaction of the neoadjuvant docetaxel drug with the DAXX and FGR1 proteins, which is consistent with the observed accumulation of cytoplasmic DAXX and overexpression of FGR1 precursors in cancer cell lines. This indicates that the DAXX and FGR1 proteins could be potential targets for the neoadjuvant docetaxel drug. The work, therefore, provides a new route for understanding the effect of the drug mode of action from the viewpoint of the change in the topology of gene-gene regulatory networks and provides a new avenue for bridging the gap between gene-gene interactions and protein-protein interactions which could have deep implications on mainstream drug development protocols.


Introduction
Recent advances in high-throughput DNA microarray technologies have fostered a growing interest in genomics, 1 proteomics, 2 and drug discovery, 3 offering a platform for a deeper understanding of gene-gene interaction and shedding light on previously uncharted avenues for understanding protein-protein and drug-protein interactions. However, a framework for understanding the causal relationship between gene-gene interactions and protein-protein interactions (PPIs) is yet to be developed. In fact, a growing number of studies have indicated that gene expression profiles could be functionally related to the protein expression levels. 4 For example, direct interactions between proteins were found to be directly linked to their gene-gene expression profiles. 5 Such interactions form a network that underlies the causal relationship between gene-gene interactions and PPIs. Learning gene regulatory networks (GRNs) from gene expression profiles could, therefore, shed light on the underlying PPI networks which could open new avenues for posterior cellular systems investigations.
Using GRNs as a route for mapping the underlying PPI networks has been the focus of a growing number of recent investigations. In a recent study, 6 it was indicated that integration between gene expressions and PPI networks could lead to prioritizing and ranking the genes most likely to be associated with breast and lung cancers. Moreover, experimental studies that complementarily use inferred GRNs and PPI networks have been able to suggest missing gene-gene interactions that does not initially show up in the inferred GRNs. 7 Currently, different approaches are being pursued to learn gene-gene regulatory networks. Many of these approaches focus on establishing methodologies for statistically unveiling correlated pairs of genes. 4,5,[8][9][10] An inherent assumption in these approaches is that if 2 genes show statistical correlation, then it is likely that they influence each other at the cellular level. In fact, coregulation relationships between genes were shown to be related to similarity in their expression profiles. 4 Machine learning of graphical models is usually used to unveil correlations between multiple genes, inferring gene-gene regulatory networks. Typical machine learning models that are widely used comprise Boolean networks, 4 Bayesian networks, 8 dynamic Bayesian networks, 9 and dependency networks. 10

Cancer Informatics
Preference of a particular graphical model over another is largely dependent on the problem under investigation. However, learning of gene-gene regulatory networks from gene expression data, that are usually sparse, often suffers from overfitting problems, leading to many pairs of genes showing correlation by chance even though they are not biologically related. 4 This severely restricts the applicability and utility of machine learning of graphical models techniques in the context of gene-gene regulatory networks.
Interestingly, incorporation of biological prior knowledge in machine learning of graphical models has shown the prospects of efficient learning of gene-gene regulatory networks such as to overcome problems incurred by sparse gene expression data. [11][12][13][14] In this context, incorporation of biological prior knowledge amounts to supplementing the graphical model technique with available information about gene expression data using cell signaling pathways relevant to the problem under investigation. This leads to restricting the variable space to lower dimensions and thereby circumventing the overfitting problems incurred in dealing with sparse gene expression profile data. Information about cell signaling pathways is available from many biological knowledge databases, such as KEGG, 15,16 CPDB, 17 REACTOME, 18 GOLD.db, 19 and PROTEIN LOUNGE. 20 However, one shortcoming of the signaling pathways available, for example, from KEGG, 15 is that they are usually represented in terms of family-to-family connections between genes rather than individual gene-to-gene connections. For example, the KEGG database shows that the MAPK kinase signaling pathway involves interaction between the fibroblast growth factor (FGF) family (22 genes) and the FGF receptor (FGFR) family (4 genes) without any reference to which individual genes are responsible for this interaction (Figure 1). To the best of our knowledge, there is no existing methodology for unveiling the specific gene-to-gene connections between gene families that underlie the generic connections represented in KEGG signaling diagrams/pathways. This definitely restricts the usability of the biological knowledge expressed in KEGG pathways for learning gene-gene regulatory networks.
Hence, the purpose of this work is 2-fold: first, to show that using family-to-family biological knowledge available from KEGG pathways along with gene expression profiles available from microarray experiments could be used to unveil individual gene-gene interactions whereof gene-gene regulatory networks could be robustly constructed and second, to use the gene-gene regulatory networks to infer the drug mode of action at the PPI level. To achieve this, we use microarray gene-gene expression profiles from 2 breast cancer cell line samples: resistant and sensitive to neoadjuvant docetaxel drug. We conduct 3 stages of investigation. First, we incorporate prior biological knowledge (family-to-family gene interactions) from the MAPK signaling pathway available from KEGG knowledge database in machine learning and feature selection graphical models to construct gene-gene regulatory networks for both samples. The reason behind choosing MAPK signaling pathway is that because we found that the gene expression profiles in our study are enriched with genes from MAPK signaling pathway (240 genes). Second, we conduct network-network comparison of the resistant and sensitive gene-gene regulatory networks to identify individual gene-gene interactions that were affected by the drug, that is, disappeared from the resistive sample. The affected gene-gene interactions are then mapped to the corresponding encoded protein chains. Third, we identify potential drug-protein interactions that underlie affected gene-gene interactions. To achieve this, we investigate the likelihood of interactions of expressed protein chains and the possibility of disruption of these interactions by the drug using protein-protein and drugprotein docking molecular modeling protocols.

Inference of gene-gene regulatory network
Affymetrix breast cancer cell line data that were subjected to the neoadjuvant docetaxel anticancer drug were retrieved from a previous study. 21 The data comprise 2 sets of samples: 14 samples that were found to be resistant to the anticancer drug and 10 samples that were found to be sensitive to it. To get meaningful gene expression data sets from the cell lines, the data (24 samples) were normalized using the Robust Multichip Average (RMA) algorithm. 6 After normalizations, the genes in the 2 data sets were annotated by their probe-IDs. The probe-IDs were obtained from the hgu95av2.db database available in Bioconductor database, which corresponds to the microarray chips used in the original study. 21 The 2 annotated data sets were then mapped to KEGG database using the KEGG.db package available in Bioconductor. This showed that many genes in the 2 data sets are annotated to MAPK signaling pathway ( Table 1). The MAPK signaling pathway was therefore used as a prior biological knowledge to Inference of gene-gene regulatory networks under Akaike information criterion-lasso restraints. We modeled the gene-gene regulatory network for each data set using Bayesian graphical models in which the graph represents the probabilistic conditional dependence between the graph vertices, genes in our case. Learning Bayesian graphical models can be achieved by assigning to each vertex (X) a number of parents (pa) and computing the corresponding conditional distributions P as follows 7 : In this context, for each gene (X), in the data sets (resistive and sensitive; Table 1), we search for a subset of (causals) parents (x i ) that best predict that gene using a linear regression model of the form: One way to find the best set of parents is to use a score function such as Akaike information criterion (AIC). 22 For a given number of parents (n), the AIC uses the residual sum of squares as the likelihood estimate that these parents are good predictors of a particular gene and incorporates a complexity penalty parameter (2p) that increases with the number of parents thereby discouraging overfitting: In fact, the AIC was shown to give more reliable results compared with alternative methods such as leave-one-out crossvalidation method (LOOCV) especially in gene expression data sets as well as in big data sets. 11 The AIC score function in that form focuses on finding best-fit parents, which does not necessarily incur concurrent best estimate of the regressive parameters (βs), namely, model selection. To alleviate this concern, the lasso estimate 23,24 was used alongside the AIC such as to ensure that the sum of the absolute values of the model regressive parameters (βs) is below a prespecified threshold parameter (s) using the following penalty function: Optimization of gene-gene regulatory networks using KEGG prior knowledge restraints and feature ranking. Systematic optimization of the value of the (s) parameter could be prohibitive due to the large dimensionality of the search space defined by the model regressive parameters (βs). To circumvent this problem, we use prior knowledge of gene-gene interactions to guide the optimization process such that all possible subsets of parents for any particular gene are restricted to only those genes that are known to interact with that gene. To achieve this, we used prior knowledge of family-to-family gene connectivity shown in the MAPK signaling pathway, available from the KEGG database. We restricted all possible subsets of parents of any particular gene to the family of genes that show connectivity with the family of that gene. For example, the MAPK signaling pathway ( Figure 1) shows that the FGF family of genes has a direct connection with the FGFR family, which indicates that the search space for predicting any of the genes in the FGFR family should be restricted only to the genes in the FGF family. A notable advantage of this procedure is that individual genegene connectivity that was not initially visible in the family-tofamily connections in the MAPK signaling pathway becomes visible. After restricting the search space for each gene using prior knowledge from MAPK signaling pathway, the univariate filtering feature selection method was used to order the search space for each gene. The filtering feature selection method is a ranking method that incrementally orders all subset of predictors according to correlation coefficients from the highest to the lowest, and therefore it is usually called feature ranking. 21 We then test for a particular value of (s) to choose the best subset of parents using the AIC score function under lasso restraint, validating the results using LOOCV. The lowest correlated parent is then iteratively removed (Algorithm 1). Using the feature ranking method with AIC score function, the lasso estimate allows, therefore, to examine different values of (s) for all possible parents of a particular gene. This procedure was found to allow for a more relaxed (s) parameter (Table 2) compared with that using only the lasso estimate described above.

GRN comparison and identifying potential drug-protein interactions
Network comparison and gene-to-protein mapping. The genegene regulatory networks for the 2 data sets, resistive and sensitive to the neoadjuvant docetaxel drug, were compared based on their pairwise gene-gene connectivity. The gene-gene

Cancer Informatics
interactions that disappeared from the resistive network were retained and are referred to as "disappeared connections" throughout. Mapping of the genes corresponding to the disappeared connections to the corresponding encoded protein chains was performed through the Human Gene Database (http://www.genecards.org/) and the corresponding structures of protein chains were obtained from the Protein Data Bank (http://www.rcsb.org/). This leads to mapping each gene to several PDB IDs that represent the same encoded protein, yet in different conformations.
Examining interaction between protein pairs corresponding to disappeared gene-gene connections: protein-protein docking. To investigate the physical reasons behind gene-gene correlation corresponding to disappeared connections, pairwise interactions of the protein chains encoded by these genes were investigated. Investigation of PPIs was performed using the ClusPro 2.0 protein-protein docking Web server (http://cluspro. bu.edu/). 25 The ClusPro server was recently reported to outperform the best human predictor groups, to select the topranked models of PPI complexes and to reliably generate high-quality structures of these complexes from the structures of separately crystallized proteins in the absence of biological information. 26 As the mode of interaction of the protein chains is largely unknown a priori, the balanced mode for computing the interaction scores was used in ClusPro. Each PPI mode is represented by ClusPro as a cluster of structure pairs that are structurally similar (root-mean-square deviation < 1.5). The number of structures of each cluster is indicative of the breadth of the free energy valley of the PPI and is used by the ClusPro server to rank different PPI modes. The scores of the top 10 interaction modes were retained.
Examining interaction of neoadjuvant docetaxel drug with protein pairs: drug-protein docking. Blind docking 27,28 of the neoadjuvant docetaxel drug was then performed against all of the protein pairs that were retained from protein-protein docking using the AutoDock 4.1 suite of programs. 29 Computation of the potential energy grids required by AutoDock at appropriate resolution is quite demanding for large proteins in terms of memory storage and CPU time; therefore, the space around each protein was partitioned using a grid that encompasses the whole protein and extends outward by 40% of its extent in each direction. A grid cell extent of 64 Å was used throughout. For each grid cell, the grid potentials were computed at 0.5-Å resolution for the docking calculations. To alleviate the possibility of missing binding sites that lie across neighboring grid cells, 2 overlapping grids were used that are shifted by 20 Å in each direction. We note that partitioning the space around the protein into grid cells that are treated independently by AutoDock serves 2 purposes: first, it circumvents the prohibitive computational cost of performing blind docking to the whole of each protein in a single run; second, because the same grid cell extent is used throughout, this setup maintains a consistent ratio of the number of docking simulations to protein size, which varies across different protein chains. In all of the docking simulations, the protein was kept rigid and fixed in space, whereas the peptide was placed at a random initial position and orientation. In the drug docking, a search in the space surrounding the protein was performed using the Lamarckian genetic algorithm (GA-LS) 30 with no restriction on the drug conformation. For every arbitrary starting position of the drug, 20 hybrid GA-LS docking runs were performed using a population size of 200, a maximum number of energy evaluations of 3 000 000, a maximum number of generations of 27 000, and 300 iterations of local search. The structures of the docked drug were stored and subsequently checked for overlap with the interface regions of each protein pairs. Drug structures that overlap with at least 3 residues of the interface of neighboring proteins, using a 3.0-Å distance threshold between the CA atoms, were selected and the drug-protein complex with the lowest binding free energy was retained in each case. Figure 2 summarizes the aforementioned methodology steps that underline the approach used in this work.

Results and Discussion
Incorporating KEGG prior knowledge and construction of gene-gene regulatory networks As described in the "Methodology" section, we incorporated prior biological knowledge (gene family-to-family interactions) from the MAPK signaling pathway that we extracted from Abbreviation: AIC, Akaike information criterion. The combination method gives a better result than using the lasso estimate only.

5
KEGG knowledgebase into the machine learning and feature selection graphical model to construct gene-gene regulatory networks from the resistive and sensitive data sets ( Table 1).
The key point is that a gene can be involved in hundreds of biological functions or pathways and it is not a priori known which pathways or gene set annotations are relevant in a given context. Incorporation of prior knowledge allowed for a statistically robust estimation of pairwise correlations between individual genes annotated in the MAPK signaling pathway. These correlations were used for the construction of the gene-gene interaction network for both sensitive and resistive samples ( Figure 3A and B).

Network comparison: resistive vs sensitive
Comparison of the resistant and sensitive gene-gene interaction networks ( Figure 3A and B) reveals that the resistant network is highly connected with distinct domains of connectivity, delineated by dashed lines in Figure 3. These domains comprise one large domain in which the RAF1 and MAP2K2 genes are pivotal for its subdomain connection, a domain in which the RAS gene is highly connected and other smaller domains that are sparsely populated. Although the resistant network shows 178 edges, only 52 edges exist in the sensitive network where the network domains, compared with the resistive network, are highly disrupted with lower degree of connectivity. Clearly, the change in the network topology embodies the effect of the docetaxel drug; however, the mechanism underlying this change is largely unknown. After all, it is the interactions between the myriad of proteins encoded by the network genes that cause gene-gene correlations and it is the drug interactions with these proteins that cause cessation of gene-gene correlations. To unveil such a mechanism, we studied these interactions at the molecular level.

Assessment of the docetaxel drug potential to disrupt the network
One way to account for such change in the gene-gene network topology is to assess the extent of disruption incurred by the drug to the interaction of the proteins encoded by the network genes. Such disruption, if any, could be expansive such as to

6
Cancer Informatics prevail the network or could be localized with a large-scale domino effect (second-order effect) throughout the network. Both scenarios could lead to network disruption; for the purpose of this work, we focused on the former. We focused on the affected network edges (pairs of genes), that is, disappeared from the resistive sample, 126 edges. We mapped each pair of genes to known pairs of protein chain structures, as outlined in the "Methodology" section and performed protein-protein docking simulations, 354 simulations, to delineate their favorable modes of interaction using the ClusPro server. 25 This leads to identification of 289 protein pairs that showed favorable modes of interactions. Favorable PPIs in ClusPro are determined by the cluster size of the docked structure clusters; this criterion has been shown to outperform human predictor groups to reliably generate high-quality structures of these complexes from the structures of separately crystallized proteins in the absence of biological information. 26 To determine whether the drug is able to disrupt these favorable proteinprotein modes of interaction, we performed blind drug-protein docking to each protein, 42 191 simulations, and inspected 3 criteria: (1) favorable drug-protein thermodynamic interaction, (2) drug interaction with the protein-protein interface region is thermodynamically more favorable than other regions on the protein surface, and (3) the drug-protein interaction takes place such that the drug overlaps at least 3 residues in the hotspot regions-residues essential for PPI 31,32 -in the protein-protein interface. These 3 criteria are tailored to target 3 competitive molecular levels of interactions in the thermodynamic domain such as to rule out the integrity of PPI in the presence of the drug. In a previous work, we used the same strategy for identifying potential lead structures that can disrupt virus assembly. 33 Imposing these criteria leads to identification of 34 PPIs, namely, 34 gene-gene interactions that are potentially disrupted by the docetaxel drug ( Figure 4A). Interestingly, most of the proteins (and thereby the corresponding genes) that satisfy the above criteria show a high degree of connectivity ( Figure  4B) are spread throughout the network ( Figure 4C). This indicates that the docetaxel drug disrupts the network by concurrent attacks that target different proteins that exist in different domains of connectivity. To further understand the cellular mechanisms that are possibly disrupted by the docetaxel drug, we need to inspect the specific gene-gene connections that are most likely affected by the drug.

Docetaxel mode of action: DAXX-mediated interaction
The drug-protein interactions that underlie the disruption of the gene-gene connections were ordered in terms of the corresponding drug-protein binding energy as shown in Figure  3A. Of these interactions, drug interaction with the DAXXgene-encoded protein (that disrupts DAXX-FAS interaction) is the most thermodynamically favorable interaction ( Figure  4A) and the drug interaction with the FGR1-encoded protein (that disrupts FGFR1-FGF13 interaction) corresponds to the highest vertex connectivity ( Figure 4B). The molecular details of theses interactions are shown in Figure 5. The molecular interactions shown in Figure 5 show faithful complementarity between interacting protein surfaces and clear overlap between the drug and the hotspots on the target proteins, providing clues for the thermodynamic stability of the PPI and highlighting the drug potential to disrupt these interactions.
In agreement with our results, highlighting the DAXX-FAS connection as a potential drug target, the DAXX protein, the protein encoded by the DAXX gene, has been reported to bind specifically to the FAS death domain. 34 Moreover, accumulation of more cytoplasmic DAXX in cancer cell lines, for example, due to interaction with a drug, is believed to participate in cellular apoptosis (programmed cell death) 35 and thereby providing routes for anticancer drug development.

7
Our results also suggest that drug FGR1 interaction has the potential to disrupt FGFR1-FGF13 interaction. Overexpression of FGFs has been reported in prostate cancer malignancies such that increased FGF signaling activates multiple signal transduction pathways, all of which play a role in prostate cancer progression. 36 This agrees with the high vertex connectivity that FGR1 shows in the resistive network ( Figure 4B) and indicates that docetaxel interaction with FGR1 probably breaks the chain of reaction throughout these signaling transduction pathways (disrupting FGR1 connections in the sensitive sample).

Conclusions
We used biological prior knowledge of family-to-family gene interactions available in the KEGG database to reveal individual gene-to-gene interaction networks that underlie the gene expression profiles of 2 cell line data sets: sensitive and resistive to neoadjuvant docetaxel breast anticancer drug. Using machine learning of graphical models to infer the topology of the gene-gene interaction networks revealed that incorporation of biological prior knowledge allows for a more robust regression model optimization compared with mainstream optimization techniques that are statistically driven by AIC penalization and lasso estimation of regression parameters.
Interestingly, comparison of the topology of the 2 networks reveals that the resistant network is highly connected with 2 large domains of connectivity: one in which the RAF1 and MAP2K2 genes form hubs of connectivity and another in which the RAS gene is highly connected. On the contrary, the sensitive network was found to be highly disrupted with a lower degree of connectivity. To unveil the physical reasons  Cancer Informatics behind such change in network topology, we investigated the interactions of the neoadjuvant docetaxel drug with the protein chains encoded by the gene-gene connectivities that disappeared from the resistive sample. We performed protein-protein docking simulations of the encoded protein chains, to delineate their thermodynamically favorable modes of interaction using the ClusPro server. 25 The potential of the neoadjuvant docetaxel drug to disrupt these thermodynamically favorable modes of interaction was further investigated by estimation of the drug-protein interaction free energy at the chain-chain interface using the AutoDock program. 29 We found that the sensitive network is likely disrupted by interaction of the neoadjuvant docetaxel drug with the DAXX and FGR1 proteins which is consistent with the observed accumulation of cytoplasmic DAXX 35 and overexpression FGR1 precursors in cancer cell lines. 36 This indicates that he DAXX and FGR1 proteins could be potential targets for the neoadjuvant docetaxel drug.
The work, therefore, provides a new route for understanding the effect of the drug mode of action from the viewpoint of the change in the topology of gene-gene regulatory networks and provides a new avenue for bridging the gap between gene-gene interactions and PPIs which could have deep implications on mainstream drug development protocols.