Drug Resistance Updates 17 (2014) 64–76

Contents lists available at ScienceDirect

Drug Resistance Updates journal homepage: www.elsevier.com/locate/drup

Prediction of dynamical drug sensitivity and resistance by module network rewiring-analysis based on transcriptional profiling Tao Zeng a , Diane Catherine Wang b , Xiangdong Wang c,∗ , Feng Xu d,∗ , Luonan Chen a,e,∗ a

Key Laboratory of Systems Biology, Shanghai Institutes for Biological Sciences, Chinese Academy of Sciences, Shanghai, China Department of Medicine, Royal Berkshire Hospital, Reading, UK Fudan University Center for Clinical Bioinformatics, Zhongshan Hospital, Fudan University Shanghai School of Medicine, China d Department of Infectious Diseases, Second Affiliated Hospital, Zhejiang University School of Medicine, Hangzhou, China e Institute of Industrial Science, University of Tokyo, Tokyo 153-8505, Japan b c

a r t i c l e

i n f o

Keywords: Drug sensitivity Resistance Systems biology Network Bioinformatics

a b s t r a c t Revealing functional reorganization or module rewiring between modules at network levels during drug treatment is important to systematically understand therapies and drug responses. The present article proposed a novel model of module network rewiring to characterize functional reorganization of a complex biological system, and described a new framework named as module network rewiring-analysis (MNR) for systematically studying dynamical drug sensitivity and resistance during drug treatment. MNR was used to investigate functional reorganization or rewiring on the module network, rather than molecular network or individual molecules. Our experiments on expression data of patients with Hepatitis C virus infection receiving Interferon therapy demonstrated that consistent module genes derived by MNR could be directly used to reveal new genotypes relevant to drug sensitivity, unlike the other differential analyses of gene expressions. Our results showed that functional connections and reconnections among consistent modules bridged by biological paths were necessary for achieving effective responses of a drug. The hierarchical structures of the temporal module network can be considered as spatio-temporal biomarkers to monitor the efficacy, efficiency, toxicity, and resistance of the therapy. Our study indicates that MNR is a useful tool to identify module biomarkers and further predict dynamical drug sensitivity and resistance, characterize complex dynamic processes for therapy response, and provide biologically systematic clues for pharmacogenomic applications. © 2014 Elsevier Ltd. All rights reserved.

1. Introduction The investigation of responses to artificial signals such as drugs or drug resistance in therapy at molecular levels is challenging, but important for the understanding of complex mechanisms during drug treatment, although there are cellular responses to distinguish external and internal challenges (Welch, 1992). The on- or off-status of drug responses is determinate and manipulatable, and used for comparative studies on drug sensitivity and resistance as a control or case condition. Therapy-responsive genes have been

∗ Corresponding authors at: Key Laboratory of Systems Biology, Shanghai Institutes for Biological Sciences, Chinese Academy of Sciences, Shanghai, China (L. Chen); Fudan University Center for Clinical Bioinformatics, Zhongshan Hospital, Fudan University Shanghai School of Medicine, China (X. Wang); Department of Infectious Diseases, Second Affiliated Hospital, Zhejiang University School of Medicine, Hangzhou, China (F. Xu). E-mail addresses: [email protected] (X. Wang), [email protected] (F. Xu), [email protected] (L. Chen). http://dx.doi.org/10.1016/j.drup.2014.08.002 1368-7646/© 2014 Elsevier Ltd. All rights reserved.

intensively investigated recently (Wang et al., 2012; Ottaviani et al., 2012; Kariko et al., 2011; Tatebe et al., 2010). It has been recognized that biological functions are generally facilitated not only by individual genes or proteins, but also by interactions between networks or modules, which are responsible for dynamical behaviors and diverse functions of living organisms. The study on dynamical drug sensitivity and efficacy to therapy can elucidate the principle of drug responses at the network level (Alsford et al., 2012; Barretina et al., 2012), based on a molecular network or structures, e.g., network motifs and modules, which can benefit the identification of efficient biomarkers in pharmacogenomic applications (Sim and Ingelman-Sundberg, 2011; Hoppe et al., 2011; Torkamani and Schork, 2012; Lussier and Li, 2012). Recent attention has been diverted from the single gene to particular functions or pathways like DNA-damage response in cancer therapy (Powell and Bindra, 2009; Lord and Ashworth, 2012). This is different from conventional studies on drug sensitivity that have mainly focused on responsive genes (Duffy et al., 2011; Yuasa et al., 2011; Chan et al., 2011; Chen et al., 2009; Cohen et al., 2011).

T. Zeng et al. / Drug Resistance Updates 17 (2014) 64–76

Fig. 1. Module Network Rewiring-analysis.

65

66

T. Zeng et al. / Drug Resistance Updates 17 (2014) 64–76

However, the mechanism by which an effective response to therapy can be achieved through gene interactions or networks of biological functions still remains unclear (Barabasi and Oltvai, 2004). The mechanism for dynamical drug sensitivity and response to therapy is considered to be rooted from functional reorganization or gene module network rewiring (Luscombe et al., 2004). Such network rewiring can be used to characterize the dynamics of drug responses and resistances, capture the functional reorganization, and extract dynamical features relevant to drug sensitivity at the functional network level (Fig. 1). Such a network-based approach can be expected to identify module biomarkers and predict dynamical drug sensitivities. A new network-based model MNR was proposed as a new computational framework in the present study, on basis of spatiotemporal transcriptional profiling data and molecular interaction data. MNR can quantitatively provide the phenotype specificity of the drug response and sensitivity differed from other types of stress responses. The drug treatment can be monitored by investigating consistent modules and re-rewiring dynamics. The present study defined a module network by a set of nodes and edges, where each consistent module of interactive genes is a node and interaction path between two modules is an edge. Such a module network can be viewed as a macro-network and a graph model for the molecular interaction system at a higher level in contrast to the traditional molecular network. Each interaction within the consistent module as a robust gene group has no rewiring or consistent during the process or various conditions, although the strength may be changed during drug treatment. While, each interaction between two consistent modules may be rewired according to therapeutic drugs, protocols, schedules, doses, or deliveries. The inter-interaction may occur at certain conditions, but disappear at other conditions. Thus, the module network and rewiring can be used to characterize dynamical drug sensitivity and resistance through the global variation and local stability of the reference molecular network in a uniform manner. Studies on the transcriptional profiling of patients with Hepatitis C virus infection (HCV) receiving Interferon therapy (Taylor et al., 2008) demonstrated that MNR can reveal changes of biological functional relationships for dynamical drug sensitivity during anti-infectious or anti-cancer therapies. Consistent modules derived by MNR can clearly distinguish two new subtypes of therapeutic responders and non-responders, while differential genes expressions failed. The rewiring interactions among consistent modules/module networks also elucidated cross-talks among biological functions, and quantitatively defined sensitivities and resistances of patients to the therapy. In addition, our work also indicates that consistent modules and rewiring networks may be associated with the improvement of drug sensitivities to classify and/or predict different dynamics of the drug sensitivity for personalized medicine.

the module network by changing interaction paths, as Step 2 in Fig. 1. Datasets used for modeling MNR include the transcriptional profiling or gene expression data (GE) organized as a tensor data and denoted as D(G,P,T). Each dg,p,t element represents the gene expression value of a gene g in gene set G for some samples p in sample set P measured at time point t from time set T. D(g,P,t) stands for the expression profile of a gene g at a time point t with all samples, D(g,p,T[s,e] ) for the expression profile of a gene g of a sample p during a time period T[s,e] which starts at time point s and ends at time point e. T can be considered as a set of condition samples in addition to time-course data. The other dataset is the bio-molecular interactome named physical interaction network (PIN) (Beyer et al., 2007). It is composed of the known protein-protein interactions and TFtarget interactions collected from different databases (KEGG (Ogata et al., 1999), ITFP (Zheng et al., 2008), Tred (Jiang et al., 2007), and iRefIndex (Razick et al., 2008)). The mathematical models are represented as a network PIN = {G,E}, where G stands for the node set of gene/protein and E for the edge set of the interaction, to measure the rewiring of a network and evaluate the conditional existence of interactions in the network. Computational analysis is used to measure the weight of each edge or the association between a pair of genes involved in an interaction. The gene pair-wise association is calculated as Pearson Correlation Coefficient (PCC). The correlation between two genes for a group of samples related to a time point is noted as PCCg1,g2,t , while the correlation between genes for a series of time points related to a sample as PCCg1,g2,p . Those two correlation values are calculated as follows:



PCCg1 ,g2 ,t =

p∈P



(dg1 ,p,t − dgP

1 ,t

(d − dgP p∈P g1 ,p,t



PCCg1 ,g2 ,p =







1

,t )

2

)(dg2 ,p,t − dgP

2 ,t



(1)

(d − dgP p∈P g2 ,p,t

2

(d − dgT1 ,p )(dg2 ,p,t t∈T g1 ,p,t

(d − dgT1 ,p ) t∈T g1 ,p,t

)

2







,t )

2

− dgT2 ,p )

(d − dgT2 ,p ) t∈T g2 ,p,t

(2) 2

P = T = d /|P| and dg,p d /|T | are means of where dg,t p∈P g,p,t t∈T g,p,t gene g expressions among a group of samples {p ∈ P} or at several time points {t ∈ T}, respectively. A consistent module is defined as a molecular group, or subnetwork, where molecular interactions within each module have no rewiring during drug treatment. It implied that intra-interactions are consistent during the process or various conditions. Gene expressions and connecting weights in the module may still vary depending on conditions. A molecular network can be re-modeled as a network of consistent modules, where each node is a consistent module and each edge is a molecular interaction path between modules.

2.2. Procedure of module network rewiring-analysis 2. MNR models 2.1. Mathematical model of dynamical drug sensitivity The present study demonstrated the dynamical drug sensitivity by two aspects according to the functional reorganization model or module network rewiring model. One aspect is to determine groups or subnetworks of cooperative genes within consistent modules across different therapeutic schedules or conditions as Step 1 in Fig. 1. Consistent modules among changing networks were selected in a common molecule-group across all therapies, to represent basic and relevant biological functions of the drug sensitivity. The second is to identify collaborations, connections, or interactions between consistent modules to construct

Main procedures for the construction of the module network rewiring model or MNR are explained in Fig. 1. The consistent network is constructed by the establishment of the molecular network in each condition based on the high throughput data and the reference molecular network (Step 0). The temporal or spatial consistent subnetworks consist of each edge between any two genes within the module without rewiring during the entire process or conditions. Consistent modules are then extracted by decomposing the consistent molecular network (Step 1). The module network for each condition is then reconstructed by connecting these modules via paths or interactions (Step 2). The functional analysis of consistent modules identifies the biological significance of genes within or between modules (Step 3), and then clinical phenotypes and genotypes of the rewiring module network are

T. Zeng et al. / Drug Resistance Updates 17 (2014) 64–76

67

Table 1 Summary of experimental datasets. Data id

Sample

Time

HCV genotype

Drug

Tissue

GSE11342 GSE7123 GSE17183

20 58 30

8 6 2

1a, 1b 1 1, 2

Peginterferon-alfa2b and ribavirin Peginterferon alpha and ribavirin IFN + Rib

PBMC PBMC Liver biopsy

mapped (Step 4). MNR is used to analyze the drug treatment, integrate clinical measurements in different schedules, and explore the potential correlation between. The relevant characteristics of MNR are extracted to be potential biomarkers for the prediction of dynamical drug sensitivities in the patient, as shown in Appendices 1 and 2.

3. Dynamical drug sensitivity characterized by MNR during antiviral therapy Three time-course microarray datasets were obtained from NCBI GEO with ID: GSE11342 (Taylor et al., 2008), GSE7123 (Taylor et al., 2007) and GSE17183 (Honda et al., 2010), and were mainly summarized in Table 1. The dataset GSE11342 was fully analyzed by MNR, while others were used for the independent validation. GSE11342 contained blood gene expression profiles of twenty patients with genotype 1 hepatitis C, untreated history, without other chronic liver diseases, and received PegIntronIM at 1.5 ␮g/kg plus ribavirin administered subcutaneously once a week for 10 weeks (Taylor et al., 2008). Those patients were sampled day 1 prior to, 3, 6, 10, 13, 27, 42, and 70 after the first injection of Interferon. According to the HCV-RNA determined by qRT-PCR with a low limit of detection of 29 IU/ml, there were 11 responders and 9 non-responders at the end of day 70. Our analyzed gene expression data contained 10,342 genes of 20 patients at 8 time points after general pre-processes, including gene symbol conversion, missing value filled with zero, duplicate gene expression average, and interpolation of temporal data. Of eight temporal networks, we found 48 consistent modules and appearance-consistent modules during therapeutic procedures, as listed in Appendix Table S1. Using the consistent modules with significant KEGG pathway enrichment, module network were re-constructed for a group of responders or non-responders at eight time points, or for each patient at six time windows, respectively. Each time window contains consecutive three time points, noted as w1 = {day -1, day 3, day 6}, w2 = {day 3, day 6, day 10}, w3 = {day 6, day 10, day 13}, w4 = {day 10, day 13, day 27}, w5 = {day 13, day 27, day 42}, and w6 = {day 27, day 42, day 70}.

4. Distinguishing genotypes of dynamical drug sensitivity by consistent modules It is questionable how new genotypes of dynamical drug sensitivity were revealed by consistent module genes (CMGs), as compared with differentially expressed genes (DEGs) in previous study (Taylor et al., 2008). A few overlapped genes could be observed between CMGs and DEGs and have significantly functional enrichment on the biological process annotation as a response to the virus, as shown in Fig. 3. Differential gene expression patterns among patients with or without effective drug responses were identified by CMGs (Fig. 2C and D) but not by DEGs (Fig. 2A and B). A consistent functional group of genes as a consistent module actually has significant biological meaning, which can help to discover unknown transcriptional features of the dynamical drug sensitivity at the biological network level. The selected modules provided the information on strong correlations with the initial degree of HCV on the basis of the known clinical data (Taylor et al., 2008). Patients were divided into groups with lower or high degrees of HCV, according to the unsupervised clustering of patients and gene expressions of the consistent modules on day 1 prior to the first injection of Interferon (Fig. 2E). There was a significant difference of HCV levels between two groups and between time points (data is not shown here). The consistent modules were used as new biomarkers to clearly divide patient groups based on subtypes, although the original research failed to show the significant relationship between differentially expressed genes and diagnostic patient HCV genotypes (Taylor et al., 2008).

5. Characterizing drug sensitivity by functional organization represented in consistent modules It is important to know the effectiveness of consistent modules in drug responses as well as responsive functions relevant to functional reorganization during dynamical drug responses. Biological meanings of consistent modules and appearance-consistent modules across therapeutic time points (Fig. 4) were evaluated by the functional pathway enrichment analysis in g:Profiler (Reimand

Table 2 Relevant cases of biological enrichment analysis of consistent modules. Module-ap-0-2

CYP2C9, CYP3A43, UGT1A10, GSTM5, UGT1A9, GPX3, UGT1A8, CYP2A6, UGT1A4, UGT1A5, UGT1A6, UGT1A7, UGT1A1, UGT1A3

Drug metabolism – other enzymes (KEGG:00983), Steroid hormone biosynthesis (KEGG:00140), Porphyrinandchlorophyll metabolism (KEGG:00860), Drug metabolism – cytochrome P450 (KEGG:00982), Metabolismofxenobiotics by cytochrome P450 (KEGG:00980), Cytochrome P450 – arranged by substrate type (REAC:211897)

Module-ap-0-1

RPS13, RPS12, RPS11, RPS10, RPS17, RPS16, RPL23A, RPS14, RPL7, RPL27A, RPS18, RPL3, RPL37A, RPS3, RPS2, RPL23, CFL1, TPT1, EIF3H, RPS15, RPL24, RPL5, RPL22, H3F3A, RPL10A, RPL4, RPS24, RPLP0, RPS23, RPS20, RPLP1, RPS29, RPL14, RPL17, RPL10, RPL12, EEF1B2, RPL19, RPL28, RPL38, RPL37, RPL34, RPL35, RPL32, RPL30, RPL31, RPS3A, UBC, EEF1G, PAK2

Influenza Infection (REAC:168254), Influenza Life Cycle (REAC:168255), Influenza Viral RNA Transcription and Replication (REAC:168273), Viral mRNA Translation (REAC:192823), Viral Protein Synthesis (REAC:192841)

Module-ap-0-18

TNFSF10, JAK2, IFNGR2

Jak-STAT signaling pathway (KEGG:04630), Influenza A (KEGG:05164)

Module-ap-1-16

CDKN1A, CCND2, STAT2

Jak-STAT signaling pathway (KEGG:04630), p53 signaling pathway (KEGG:04115), Cell cycle (KEGG:04110), Hepatitis C (KEGG:05160), Transcriptional activation of p53 responsive genes (REAC:69560)

68

T. Zeng et al. / Drug Resistance Updates 17 (2014) 64–76

Fig. 2. Differences between consistent module genes and differentially expressed genes in Data GSE11342. The expression profiles of differentially expressed genes for responders (A) or non-responders (B) on eight days can be compared with those of consistent module genes for responders (C) or non-responders (D). The clustering of patients based on consistent module profiles on day 1 prior to the first injection of Interferon reveals two possible subtypes of patients corresponding to different clinical phenotypes, especially related to HCV initial degree (E).

T. Zeng et al. / Drug Resistance Updates 17 (2014) 64–76

69

Fig. 3. Differences between consistent module genes and differentially expressed genes in Data GSE11342. The known biological networks with consistent module genes are compared with differentially expressed genes. Yellow nodes represent for consistent module genes (CMGs), green nodes for differentially expressed genes (DEGs), and red nodes for the overlap.

Fig. 4. Network of consistent modules and appearance-consistent modules across the whole therapy period. (A) The network of consistent modules. (B) The network of appearance-consistent modules.

70

T. Zeng et al. / Drug Resistance Updates 17 (2014) 64–76

et al., 2007). The significant modules (p < 0.05) were used to build static module networks and investigate dynamical characteristics. The pathways enriched in consistent modules suggested that several pathways play consistent roles in dynamical drug responses during the therapy, as shown in Table 2. The complete table of modules and functional enrichments was listed in Appendix Table S1. For example, pathways of drug metabolism-enzymes (KEGG:00983) and cytochrome P450 (KEGG:00982) were observed in module Module-ap-0-2, while pathways of Viral mRNA Translation (REAC:192823) and Viral Protein Synthesis (REAC:192841) were enriched in module Module-ap-0-1. It indicates that the HCV therapy response may have two global characteristics as the inputs and outputs of the drug therapy, where viral genetic events can be influenced by the drug therapy in a drug-specific metabolic way. The drug sensitivity of each patient can be determined according to the consistent modules and characteristics. The relationships between a signaling pathway and a virusrelevant pathway in Interferon therapy (Taylor et al., 2008) were also extracted, such as Influenza A (KEGG:05164) and JakSTAT signaling pathway (KEGG:04630) in module Module-ap-0-18, and Hepatitis C (KEGG:05160) and Jak-STAT signaling pathway (KEGG:04630) in module Module-ap-1-16. It suggests that virusused pathways and general signaling pathways may share key genes in multiple diseases. The nature of drug sensitivity during the therapy could be associated with the regulation of the signaling pathway to control relevant genes (Taylor et al., 2008), in order to prevent the potential activation of disease-specific pathways. The cross-talking between signaling pathways may also be responsible for drug sensitivity because the gene overlap among Hepatitis C (KEGG:05160), Jak-STAT signaling pathway (KEGG:04630), and p53 signaling pathway (KEGG:04115) also appears in module Module-ap-1-16. The treatment with interferon for HCV is a type of gene therapy through the Jak-STAT signaling pathway (Taylor et al., 2008), while the challenge to be faced during the application is the p53 toxicity to normal tissues caused by the pharmacological p53 activation (Vassilev, 2007). It indicates that the therapy per se can cure some kinds of diseases, but may simultaneously cause side-effects through complex and unknown cross-talks within toxic pathways. Consistent modules might be involved in different biological mechanisms of dynamical drug sensitivity and resistance, including drug metabolism pathways, drug-target signaling pathways, and complicated pathway crosstalks. Therefore, the dynamical expression states of those modules and interactions play critical and dominant roles in the regulation and control of biological networks.

6. Revealing global functional relationship in effective responses to therapy by quantitative MNR The module genes keep consistent connections and interactions during the therapy, and contribute to several disease- or therapyassociated pathways. It is worth to further investigate the activity variance of modules and rewiring connections during the drug therapy, which can disclose differences of the drug sensitivity and resistance to therapy even among patients. The average of differential expression E-score (fold-change here) (Taylor et al., 2008) and newly defined expression-correlation scores EC-score (as module activity score defined in formula (6) in Appendix) indicate that a few modules marked with star in Fig. 5A and B had significant differences between responders and non-responders. Interaction weights between consistent module pairs as a module interaction activity score defined in formula (7) in Appendix could distinguish between patients with or without responses to therapy, based on the most discriminative consistent module pairs and interactions shown in Fig. 5C. Elements from patients with or

Fig. 5. E-score/EC-score of consistent modules and connection weights of consistent module pairs. The E-score profiles (A), EC-score profiles (B), weight profiles (C) of consistent modules were recorded in patients with or without responses for seven therapy days, of which the values on the day one prior to the first injection of Interferon were used as control. The heat-maps demonstrate that each column represents a group of patients, where R-ID and NR-ID stand for patients with or without responses, respectively. Each row represents a module interaction with the label as Mx –My , where Mx and My are custom ID of two modules involved in an interaction.

T. Zeng et al. / Drug Resistance Updates 17 (2014) 64–76

without therapeutic responses tended to construct close module interactions during the drug therapy, while responders had the stronger module connections (Fig. 5C). It indicates that functional modules need to communicate at certain time periods of therapy. The drug sensitivity may be non-effective for unresponsive patients due to the lack of necessary genes to shorten and/or bridge connections of modules on the biological network. The interactions between modules and the temporal rewiring may reflect the close responsive relationships among disease-specific and drug-specific functions enriched in modules. Recent studies demonstrated that the global characteristic of signaling pathways represented by network entropy was correlated with the survivability of patients with or without cancer (Breitkreutz et al., 2012). However, there are a number of limitations when the single pathway was studied in a static manner. It also failed to reflect the variance of the association in a dynamical manner. The MNR and related module interaction weights represented the sum-weight of shortest paths in the original molecular network, and uncovered a deeper association between larger biological networks involved in multiple pathways and patient phenotypes of drug sensitivities and resistances. 7. Predicting responder and non-responder with drug sensitivity by spatio-temporal markers The analysis of the MNR demonstrated significant differences of biological networks between patients with or without responses.

71

The question is whether those differences determined by MNR can be used to distinguish patients. The MNR chain could start and end at any pair of check time points in a patient, after the time period was optimized to classify patients. In the experiment, the time window within consecutive three time points in the original data was used. The distinct module connection weights of patients with or without responses at certain time periods were analyzed by PermutMatrix (Caraux and Pinloche, 2005), including the early time windows (Fig. 6) to reveal systematic clues for earlier diagnosis. Module interactions were clustered into six main groups corresponding to six time windows (w1-w6), to reflect spatial (patients) and temporal (time windows) specificities simultaneously. Patients with responses had strong correlations and interactions among modules at middle time windows (w3-w5) covering days 06, 10, 13, 27 and 42 after drug injection. Clinical studies suggested that early changes during Interferon therapy be important to predict the therapeutic outcome (Farci et al., 2002). HCV patients with an early virologic response should have a better chance of being cured (Sarasin-Filipowicz et al., 2008). It implies that functional modules should have strongly correlated networks and communications for the dynamical drug sensitivity at certain time periods of the therapy. Module interactions with time windows shown in each row of the heat-map in Fig. 6 were used as features to build different classifiers (Hall et al., 2009) and the leave-one-out cross-validation was used to evaluate the classification model performance, based on

Fig. 6. Hierarchical clustering of module network rewiring chains. The horizontal and vertical dendrograms represent for the hierarchical tree of responders and nonresponders and of discriminative consistent module interactions within given time windows. Each column represents a patient, where R-ID and NR-ID stand for patients with or without responses. Each row represents a consistent module interaction in a time window with the label as Mx –My –wid , where Mx and My are custom IDs of two modules involved in an interaction and wid is the ID of time window. (A) The discriminative consistent module interactions in early time period (e.g. during the first 3 time windows). (B) The discriminative consistent module interactions in later time period (e.g. during the last 3 time windows).

72

T. Zeng et al. / Drug Resistance Updates 17 (2014) 64–76

the chain encoding vectors (i.e. numeric vector of module interaction weight in each time window). The accuracy and AUC of the rule-based prediction model as the decision tree and the non-rulebased prediction model as support vector machine or logistic model were 0.85 and 1.0, respectively. The additional results from the independent validation in the next subsection also strongly support the effect of those module biomarkers in biomedical study and application. 8. Validating the accuracy and robustness of module biomarker on independent groups of HCV patients Consistent modules or quantified module interaction weights were proposed to have the close association with drug sensitivity or to classify and predict dynamical drug sensitivities or therapy responses of patients, based on results from Data GSE11342 obtained through MNR. Furthermore, we focused on the independent validation of those two hypotheses on Data GSE7123 (Taylor et al., 2007) and Data GSE17183 (Honda et al., 2010) by consistent modules mined from Data GSE11342. Those known consistent modules were used to re-construct module networks in

independent datasets. Data GSE7123 on gene expression data of 58 patients at one time point before drug injection and five points during therapy were used to validate the association between consistent modules and drug sensitivities, and the sensitivity prediction by quantified module interaction weights; while Data GSE17183 only had gene expression data of 30 patients at time points before and after drug injection for the validation of the association between consistent modules and drug sensitivity. Taylor et al. (2007) and Huang et al. (2008) analyzed patient samples in Data GSE7123, including African-American patients, and Caucasian-American patients, respectively. Honda et al. (2010) analyzed the samples of liver biopsy, liver lobules, or cells in portal areas in Data GSE17183, respectively. We made three replicated validations of different populations and tissues in Data GSE7123 and Data GSE17183, respectively. Fig. 7 demonstrated the value distribution of weights of possible module interactions from the analysis of box-and-whisker diagram in patients with or without response at a time point in one of three datasets, including the distribution of weights of responders and non-responders corresponding to the same patients at eight time points from Data GSE11342 (Fig. 7A), to three patient populations at six time points

Fig. 7. Distribution of interaction weights between consistent modules for responders and non-responders in multiple datasets. It demonstrates the weight distribution of Data GSE11342 (A), Data GSE7123 corresponding to patients (ALL, AA, CA) (B)–(D), or Data GSE17183 to three different patient tissues (LB, CLL, CPA) (E)–(G). R.ID and NR.ID represent patients with or without responses.

T. Zeng et al. / Drug Resistance Updates 17 (2014) 64–76 Table 3 Prediction performance of therapy response on Data GSE7123. Method Huang et al. MNR

All patients 72.4% 81%

AA patients 85.7% 100%

73

This work was also supported by grants from Natural Science Foundation of Zhejiang Province (LR12H01003). CA patients 100% 90%

Appendix A. A.1. Module Network Rewiring-analysis

from Data GSE7123 (Fig. 7B–D), to three tissues at two time points from Data GSE17183 (Fig. 7E–G), respectively. Patients with responses had less module interaction weights than those without responses, and significant differences between two numeric vectors of module interaction weights or functional associations at days during drug treatment (Fig. 7A). MNR-based biomarkers were used to build the classification model with the decision tree, SVM, or logistic model, which had higher prediction accuracy as compared to previous methods, as shown in Table 3. The performance was achieved by consistent modules mined in Data GSE11342 and validated on Data GSE7123, which implies that such consistent modules are robust and accurate in the understanding of dynamical drug sensitivities in patients receiving Interferon therapy for hepatitis C treatment. Those consistent modules relevant to anti-infection therapy were also validated in cancer samples (e.g. HCV-related Hepatocellular Carcinoma) and would be helpful in designing anti-cancer therapy, seen in Appendix 3. 9. Conclusion The proposed MNR from the present study analyzed rewired associations between consistent functions based on the module network rewiring, different from molecular interactions or dysfunctions with the rewired interactions of molecules. MNR can analyze the dynamic drug sensitivity or resistance of the functional interaction rather than the molecular interaction. Several consistent modules and varying associations were found to distinguish between patients with or without drug responses by MNR on the gene expression data from HCV patients receiving Interferon injection. MNR can also be applied to identify edge biomarkers (Yu et al., 2013; Zhang et al., 2014), network biomarkers (He et al., 2012; Liu et al., 2012; Zeng et al., 2013), and even dynamical network biomarkers (Chen et al., 2012; Li et al., 2013; Liu et al., 2013; Zeng et al., 2014), and track dynamics of drug responses or resistances of anti-infection/anti-cancer therapies, or even phenotype evolutions of general dynamical biological processes. Our study indicates that MNR is a powerful tool to quantitatively predict dynamical drug sensitivities and resistances by identifying consistent modules, characterize complex dynamic processes for therapeutic responses, and provide biologically systematic clues for pharmacogenomic applications. Acknowledgements This work was supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (XDB13040700), and the National Program on Key Basic Research Project [Grant No. 2014CB910504]. This work was supported by National Natural Science Foundation of China (NSFC) with Nos. 61134013, 91029301, 31200987, 91130033, 91230204, 81270099, 8132018001, 81270131, and 81370176. This work was also supported by Knowledge Innovation Program of SIBS, CAS with Grant No. 2013KIP218. The work was supported by Zhongshan Distinguished Professor Grant (XDW), Shanghai Committee of Science and Technology (12JC1402200, 12431900207, 11410708600), Zhejiang Provincial Natural Science Foundation (Z2080988), Zhejiang Provincial Science Technology Department Foundation (2010C14011), and Ministry of Education, Academic Special Science and Research Foundation for PhD Education (20130071110043).

Actually, many previous studies have discussed the cooperative genes as modules or the varying gene cooperation as gene network rewiring, where a module is generally considered as a basic unit to perform biological functions at the network level. Different from mining the conserved gene expression (Herschkowitz et al., 2007), discovering modules is a difficult computation task, and many methods have been applied to this problem (Langfelder et al., 2011; Hou et al., 2011; Komurov and White, 2007; Azuaje et al., 2010; Beltrao et al., 2010; Jaimovich et al., 2010), but there is few work to detect the consistent modules across different conditions, and also few researchers have exploited dynamic interactions among modules to analyze network functions although there are some simple schemes using the overlapped genes to link different modules (Ihmels et al., 2002; Liu et al., 2010; Kovacs et al., 2010). In addition to detecting modules, study on gene network rewiring behind cellular response has also attracted much attention and becomes a key biological model in recent studies (Shou et al., 2011; Sun et al., 2012; Ihmels et al., 2005; De Smet and Van de Peer, 2012; Freschi et al., 2011; Roguev et al., 2008; Bandyopadhyay et al., 2010; Bhardwaj et al., 2010). However, these methods mainly paid attention on characterizing the changes between molecules, rather than dynamical regulations of functional modules during spatio-temporal rewiring processes. In order to reveal the essential mechanism of functional reorganization at the network level, it is necessary to develop a theoretical model to combine the module and network rewiring together in a dynamic and unified manner. A.1.1. Construction of biomolecular network (Step 0 in Fig. 1) The molecular networks are constructed based on the available data (see Table 1, Step 0 in Fig. 1, and Result and Discussion), e.g., expression data and the reference/background network. There are two kinds of context-based networks focused in this study. One is the time-specific network (across patients or samples) at one time point or a time window, and the other one is the patientspecific (sample-specific) network (across whole time periods) for one patient or a patient group. In this step of network construction, the correlations between two genes are calculated. For the patientspecific network, PCCg1,g2,p is used, which is based on patient’s expression values during a time window (i.e., several time points in the time window, or values of several replicates at a given time point). While, for the time-specific network, PCCg1,g2,t is adopted, which is based on a group of patients’ expression values at a given time point (Note that the time-specific network for one sample will be trivial and equivalent to a patient-specific network at the same time point ideally). According to the critical values of correlations (Vaughan, 2001), the edges in PIN (protein interaction network) having correlations less than the statistic thresholds (or P-value >0.05) are removed and the remaining edges consist of the so-called PIN co-expression networks (Zeng and Chen, 2012) modeling the context-based networks. For convenience, the patient-specific network is denoted as PIN(p,T), where p is a patient and T is a time period. And the time-specific network is denoted as PIN(P,t), where t is a time point and P is a group of patients. A.1.2. Decomposition of biomolecular network by consistent modules (Step 1 in Fig. 1) The molecular network is decomposed based on edge consistency across samples and times. In order to find consistent

74

T. Zeng et al. / Drug Resistance Updates 17 (2014) 64–76

subnetworks across multiple context-based networks (e.g. timespecific networks here), a consecutive version of frequent graph mining method (Han et al., 2007) is adopted. The overlapped edges presenting in the networks {PIN(P,t)}t∈T consist of the socalled consistent network S. Similarly, the overlapped edges only appearing in the networks after particular time points consist of the appearance-consistent network Sk , while the overlapped edges only existing in the networks before particular time points consist of the absence-consistent network Sk . Formally, the three kinds of consistent networks can be obtained according to following definitions. S=



t∈T

Sk = Sk =

 

(3)

PIN(P, t) PIN(P, t) −

t∈T, t≥k

PIN(P, t) −

t∈T, t 0|x ∈ M , y ∈ M |}

 (M , M ) ∈ E(tp )

(7)

where ex is the differential expression of gene e, and sp(x,y) is the length of the shortest route between genes x and y on the context-based molecular network, and sp(x,y) equals to the given maximum value when x is not connected to y. Note that WV is related to the correlation but WE is related to the distance between the two components. A module is considered to be more active or more functional in the biological system when

its members have close correlation and high differential expression simultaneously so that the weight WV is large. Otherwise, the module is inactive. On the other hand, the weight of modules’ association reflects the closeness among all genes in the two modules. A small weight WE (distance) between two modules means that the two modules may be the two components of another function module on a high level of the biological network because of their strong functional correlation. Otherwise, the two modules will be independent when the weight of their connection is large. A.1.4. Functional annotation of consistent modules (Step 3 in Fig. 1) The module network rewiring aims to study the functional reorganization, rather than conventional molecule network rewiring. For each module, its functions are re-analyzed by conventional functional enrichment analysis (Subramanian et al., 2005) or network ontology analysis (Wang et al., 2011). A.1.5. Spatial and temporal module network rewiring (Step 4 in Fig. 1) At one time window (as consecutive three time points in data here), different samples/patients have respective module networks with the same nodes (modules) but different edges (module connections), which are resulted from the spatial module network rewiring and can reflect the patient specificity. Meanwhile, for one patient, there is also a sequence of module networks with the same nodes but different edges among them at different time windows. These changes of the network structures in a temporal  order  are p p defined as the module network rewiring chain CT = t for t∈T the corresponding patient, which can describe the dynamics of a module network (e.g. temporal module network rewiring) during a biological process like therapy response and potentially reflect the patient’s response genotype. In fact, the spatial module network rewiring can be considered as the trivial condition of the difference  p p among multiple module network rewiring chains (C{k} = t ) t≡k when the chains’ length is equal to one. Therefore, the module network rewiring chain can be used to unify the spatial and temporal module network rewiring. A.2. Experiment of module network rewiring-analysis (1) The temporal consistent networks without the patientspecificity are identified (based on all patients’ gene expression at each time point) and the consistent modules are naturally detected from maximal network components, under the assumption that a therapy is mainly to affect the association among functions represented by modules. (2) Then, for the two patient groups with therapy response or not, their module networks are re-built again at different time points in order to trace the network rewiring with the patientresponse (patient group) specificity. This is a static strategy by using the snapshot of the module network to distinguish patient groups with different responses. (3) For a dynamical strategy, every patient can also have several structures of module networks corresponding to different time windows (covering consecutive three time points) during the therapy period. The module network rewiring chain for any patient is expected to reflect patients’ responses to the therapy in a dynamic manner, so that each chain is encoded to a (WE ) weighted feature vector similar to 0-1 vector shown in Fig. 1. (4) Finally, traditional hierarchical clustering is used to show the categories of the module network rewiring chains corresponding to patients, and common machine learning methods are

T. Zeng et al. / Drug Resistance Updates 17 (2014) 64–76

used to build appropriate prediction model for judging drug therapy effect on patients at early time point. A.3. Exploring the active role of consistent modules in reverse translational research To be a main merit in translational research, our findings can provide novel network based module biomarker to predict dynamical drug sensitivity, which can help physicians to select an optimal treatment for individual patient. Meanwhile, these markers also have ability to assist fundamental biological mechanism research like disease development and progression focused in reverse translational research (Sinha et al., 2011). HCV is a well-known relevant factor to cause the deterioration of liver disease as liver cancer. The gene relationships between HCV and HCV-induced HCC cohorts are attracting more attention. Thus, the fourth dataset from GEO with ID as GSE17967 is used to validate the extensive application of consistent modules, which contains gene expressions of PBMC samples from HCV cirrhotic tissues with (16 HCV-HCC) and without (47 HCV) HCC (Archer et al., 2009). For our mined consistent modules above, the numbers of significantly recovered modules in HCV cohort (context-based network for HCV samples) are 27 (75% recovered ratio means a module has 75% interactions recovered), 39 (50% recovered ratio), and 43 (25% recovered ratio). While, the numbers of significantly recovered modules in HCV-HCC cohort (context-based network for HCV-HCC samples) are 9 (75% recovered ratio), 22 (50% recovered ratio), and 28 (25% recovered ratio). Therefore, our consistent modules show their robustness on HCV but not HCV-HCC relevance again in this independent data. More importantly, there are 4 consistent modules are still completely observed in HCV-HCC cohort: Module-ap-1-18 (PSEN1, PAK1, CASP1); Module-ap-0-27 (HLA-C, HLA-B, HLA-A); Moduleap-1-15 (SCP2, ACOX1, HSD17B4); and Module-ap-1-8 (H2BFS, HIST1H2BF, HIST2H2AA3, HIST1H2BK). Although few of them have significant enrichment on known biological pathways, they are indeed related to HCC according to the HCCnet (He et al., 2010). This expert database (April 18, 2010) contains 2234 HCC-related genes collected from multiple datasets and supplies many summaries on each HCC-related genes, which include all the genes of 4 consistent modules. Table S2 lists each gene’s two statistics from HCCnet: the gene ontology and expression alteration in each fundamental dataset. It is easy to see that these genes are actually different from DEGs with extremely up-regulation or down-regulation, and they tend to have varying expression alteration (but potential consistent expression associations). As dys-up-regulated genes are always interested by biologist or clinician, it is noticeable that PAK1 in module Module-ap-1-18, and all four genes H2BFS, HIST1H2BF, HIST2H2AA3, HIST1H2BK in Module-ap-1-8 are consistently up-regulated in different HCV induced hepatocellular carcinoma or HCV-associated hepatocellular carcinoma samples (Table S2). These dys-up-regulated consistent module genes indeed have found respective roles in HCV or HCC (Ching et al., 2007; Ishida et al., 2007; Lu et al., 2007; Miura et al., 2008). Besides, all the consistent module genes also show better classification performance than DEGs (CMGs and DEGs are both from GSE11342) on the hierarchical clustering of HCV and HCV-HCC samples (from GSE17967) (seeing Figure S1). Therefore, our findings provide new candidate targets in human pathogen study, and prepares next round of translational research well. Appendix B. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/ j.drup.2014.08.002.

75

References Alsford, S., Eckert, S., Baker, N., Glover, L., Sanchez-Flores, A., Leung, K.F., Turner, D.J., Field, M.C., Berriman, M., Horn, D., 2012. High-throughput decoding of antitrypanosomal drug efficacy and resistance. Nature 482 (7384), 232–236. Archer, K.J., Mas, V.R., David, K., Maluf, D.G., Bornstein, K., Fisher, R.A., 2009. Identifying genes for establishing a multigenic test for hepatocellular carcinoma surveillance in hepatitis C virus-positive cirrhotic patients. Cancer Epidemiol. Biomarkers Prev. 18 (11), 2929–2932. Azuaje, F., Devaux, Y., Wagner, D.R., 2010. Coordinated modular functionality and prognostic potential of a heart failure biomarker-driven interaction network. BMC Syst. Biol. 4, 60. Bandyopadhyay, S., Mehta, M., Kuo, D., Sung, M.K., Chuang, R., Jaehnig, E.J., Bodenmiller, B., Licon, K., Copeland, W., Shales, M., et al., 2010. Rewiring of genetic networks in response to DNA damage. Science 330 (6009), 1385–1389. Barabasi, A.L., Oltvai, Z.N., 2004. Network biology: understanding the cell’s functional organization. Nat. Rev. Genet. 5 (2), 101–113. Barretina, J., Caponigro, G., Stransky, N., Venkatesan, K., Margolin, A.A., Kim, S., Wilson, C.J., Lehar, J., Kryukov, G.V., Sonkin, D., et al., 2012. The Cancer Cell Line Encyclopedia enables predictive modelling of anticancer drug sensitivity. Nature 483 (7391), 603–607. Beltrao, P., Cagney, G., Krogan, N.J., 2010. Quantitative genetic interactions reveal biological modularity. Cell 141 (5), 739–745. Beyer, A., Bandyopadhyay, S., Ideker, T., 2007. Integrating physical and genetic maps: from genomes to interaction networks. Nat. Rev. Genet. 8 (9), 699–710. Bhardwaj, N., Kim, P.M., Gerstein, M.B., 2010. Rewiring of transcriptional regulatory networks: hierarchy, rather than connectivity, better reflects the importance of regulators. Sci. Signal 3 (146), ra79. Breitkreutz, D., Hlatky, L., Rietman, E., Tuszynski, J.A., 2012. Molecular signaling network complexity is correlated with cancer patient survivability. Proc. Natl. Acad. Sci. U. S. A. 109 (23), 9209–9212. Caraux, G., Pinloche, S., 2005. PermutMatrix: a graphical environment to arrange gene expression profiles in optimal linear order. Bioinformatics 21 (7), 1280–1281. Chan, E., Prado, D.E., Weidhaas, J.B., 2011. Cancer microRNAs: from subtype profiling to predictors of response to therapy. Trends Mol. Med. 17 (5), 235–243. Chen, B.J., Causton, H.C., Mancenido, D., Goddard, N.L., Perlstein, E.O., Pe’er, D., 2009. Harnessing gene expression to identify the genetic basis of drug resistance. Mol. Syst. Biol. 5, 310. Chen, L., Liu, R., Liu, Z.P., Li, M., Aihara, K., 2012. Detecting early-warning signals for sudden deterioration of complex diseases by dynamical network biomarkers. Sci. Rep. 2, 342. Ching, Y.P., Leong, V.Y., Lee, M.F., Xu, H.T., Jin, D.Y., Ng, I.O., 2007. P21-activated protein kinase is overexpressed in hepatocellular carcinoma and enhances cancer metastasis involving c-Jun NH2-terminal kinase activation and paxillin phosphorylation. Cancer Res. 67 (8), 3601–3608. Cohen, A.L., Soldi, R., Zhang, H., Gustafson, A.M., Wilcox, R., Welm, B.E., Chang, J.T., Johnson, E., Spira, A., Jeffrey, S.S., et al., 2011. A pharmacogenomic method for individualized prediction of drug sensitivity. Mol. Syst. Biol. 7, 513. De Smet, R., Van de Peer, Y., 2012. Redundancy and rewiring of genetic networks following genome-wide duplication events. Curr. Opin. Plant Biol. 15 (2), 168–176. Duffy, M.J., O’Donovan, N., Crown, J., 2011. Use of molecular markers for predicting therapy response in cancer patients. Cancer Treat. Rev. 37 (2), 151–159. Farci, P., Strazzera, R., Alter, H.J., Farci, S., Degioannis, D., Coiana, A., Peddis, G., Usai, F., Serra, G., Chessa, L., et al., 2002. Early changes in hepatitis C viral quasispecies during interferon therapy predict the therapeutic outcome. Proc. Natl. Acad. Sci. U. S. A. 99 (5), 3081–3086. Freschi, L., Courcelles, M., Thibault, P., Michnick, S.W., Landry, C.R., 2011. Phosphorylation network rewiring by gene duplication. Mol. Syst. Biol. 7, 504. Hall, M., Frank, E., Holmes, G., Pfahringer, B., Reutemann, P., Witten, I.H., 2009. The WEKA data mining software: an update. SIGKDD Explor Newsl. 11 (1), 10–18. Han, J., Cheng, H., Xin, D., Yan, X., 2007. Frequent pattern mining: current status and future directions. Data Min Knowl. Discov. 15 (1), 55–86. He, B., Qiu, X., Li, P., Wang, L., Lv, Q., Shi, T., 2010. HCCNet: an integrated network database of hepatocellular carcinoma. Cell Res. 20 (6), 732–734. He, D., Liu, Z.P., Honda, M., Kaneko, S., Chen, L., 2012. Coexpression network analysis in chronic hepatitis B and C hepatic lesions reveals distinct patterns of disease progression to hepatocellular carcinoma. J. Mol. Cell Biol. 4 (3), 140–152. Herschkowitz, J.I., Simin, K., Weigman, V.J., Mikaelian, I., Usary, J., Hu, Z., Rasmussen, K.E., Jones, L.P., Assefnia, S., Chandrasekharan, S., et al., 2007. Identification of conserved gene expression features between murine mammary carcinoma models and human breast tumors. Genome Biol. 8 (5), R76. Honda, M., Nakamura, M., Tateno, M., Sakai, A., Shimakami, T., Shirasaki, T., Yamashita, T., Arai, K., Sakai, Y., Kaneko, S., 2010. Differential interferon signaling in liver lobule and portal area cells under treatment for chronic hepatitis C. J. Hepatol. 53 (5), 817–826. Hoppe, R., Brauch, H., Kroetz, D.L., Esteller, M., 2011. Exploiting the complexity of the genome and transcriptome using pharmacogenomics towards personalized medicine. Genome Biol. 12 (1), 301. Hou, L., Wang, L., Qian, M., Li, D., Tang, C., Zhu, Y., Deng, M., Li, F., 2011. Modular analysis of the probabilistic genetic interaction network. Bioinformatics 27 (6), 853–859. Huang, T., Tu, K., Shyr, Y., Wei, C.C., Xie, L., Li, Y.X., 2008. The prediction of interferon treatment effects based on time series microarray gene expression profiles. J. Transl. Med. 6, 44.

76

T. Zeng et al. / Drug Resistance Updates 17 (2014) 64–76

Ihmels, J., Friedlander, G., Bergmann, S., Sarig, O., Ziv, Y., Barkai, N., 2002. Revealing modular organization in the yeast transcriptional network. Nat. Genet. 31 (4), 370–377. Ihmels, J., Bergmann, S., Gerami-Nejad, M., Yanai, I., McClellan, M., Berman, J., Barkai, N., 2005. Rewiring of the yeast transcriptional network through the evolution of motif usage. Science 309 (5736), 938–940. Ishida, H., Li, K., Yi, M., Lemon, S.M., 2007. p21-activated kinase 1 is activated through the mammalian target of rapamycin/p70 S6 kinase pathway and regulates the replication of hepatitis C virus in human hepatoma cells. J. Biol. Chem. 282 (16), 11836–11848. Jaimovich, A., Rinott, R., Schuldiner, M., Margalit, H., Friedman, N., 2010. Modularity and directionality in genetic interaction maps. Bioinformatics 26 (12), i228–i236. Jiang, C., Xuan, Z., Zhao, F., Zhang, M.Q., 2007. TRED: a transcriptional regulatory element database, new entries and other development. Nucleic Acids Res. 35 (Database issue), D137–D140. Kariko, K., Muramatsu, H., Ludwig, J., Weissman, D., 2011. Generating the optimal mRNA for therapy: HPLC purification eliminates immune activation and improves translation of nucleoside-modified, protein-encoding mRNA. Nucleic Acids Res. 39 (21), e142. Komurov, K., White, M., 2007. Revealing static and dynamic modular architecture of the eukaryotic protein interaction network. Mol. Syst. Biol. 3, 110. Kovacs, I.A., Palotai, R., Szalay, M.S., Csermely, P., 2010. Community landscapes: an integrative approach to determine overlapping network module hierarchy, identify key nodes and predict network dynamics. PLoS One 5 (9). Langfelder, P., Luo, R., Oldham, M.C., Horvath, S., 2011. Is my network module preserved and reproducible? PLoS Comput. Biol. 7 (1), e1001057. Li, M., Zeng, T., Liu, R., Chen, L., 2013. Detecting tissue-specific early warning signals for complex diseases based on dynamical network biomarkers: study of type 2 diabetes by cross-tissue analysis. Brief. Bioinform., http://dx.doi.org/10.1093/bib/bbt1027. Liu, C., Li, J., Zhao, Y., 2010. Exploring hierarchical and overlapping modular structure in the yeast protein interaction network. BMC Genomics 11 (Suppl4), S17. Liu, X., Liu, Z.P., Zhao, X.M., Chen, L., 2012. Identifying disease genes and module biomarkers by differential interactions. J. Am. Med. Inform. Assoc. 19 (2), 241–248. Liu, R., Wang, X., Aihara, K., Chen, L., 2013. Early diagnosis of complex diseases by molecular biomarkers, network biomarkers, and dynamical network biomarkers. Med. Res. Rev., http://dx.doi.org/10.1002/med.21293. Lord, C.J., Ashworth, A., 2012. The DNA damage response and cancer therapy. Nature 481 (7381), 287–294. Lu, Y.S., Kashida, Y., Kulp, S.K., Wang, Y.C., Wang, D., Hung, J.H., Tang, M., Lin, Z.Z., Chen, T.J., Cheng, A.L., et al., 2007. Efficacy of a novel histone deacetylase inhibitor in murine models of hepatocellular carcinoma. Hepatology 46 (4), 1119–1130. Luscombe, N.M., Babu, M.M., Yu, H., Snyder, M., Teichmann, S.A., Gerstein, M., 2004. Genomic analysis of regulatory network dynamics reveals large topological changes. Nature 431 (7006), 308–312. Lussier, Y.A., Li, H., 2012. The rise of translational bioinformatics. Genome Biol. 13 (8), 319. Miura, K., Taura, K., Kodama, Y., Schnabl, B., Brenner, D.A., 2008. Hepatitis C virusinduced oxidative stress suppresses hepcidin expression through increased histone deacetylase activity. Hepatology 48 (5), 1420–1429. Ogata, H., Goto, S., Sato, K., Fujibuchi, W., Bono, H., Kanehisa, M., 1999. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 27 (1), 29–34. Ottaviani, D., Lever, E., Mao, S., Christova, R., Ogunkolade, B.W., Jones, T.A., Szary, J., Aarum, J., Mumin, M.A., Pieri, C.A., et al., 2012. CTCF binds to sites in the major histocompatibility complex that are rapidly reconfigured in response to interferon-gamma. Nucleic Acids Res. 40 (12), 5262–5270. Powell, S.N., Bindra, R.S., 2009. Targeting the DNA damage response for cancer therapy. DNA Repair (Amst) 8 (9), 1153–1165. Razick, S., Magklaras, G., Donaldson, I.M., 2008. iRefIndex: a consolidated protein interaction database with provenance. BMC Bioinform. 9, 405. Reimand, J., Kull, M., Peterson, H., Hansen, J., Vilo, J., 2007. g:Profiler–a webbased toolset for functional profiling of gene lists from large-scale experiments. Nucleic Acids Res. 35 (Web Server issue), W193–W200. Roguev, A., Bandyopadhyay, S., Zofall, M., Zhang, K., Fischer, T., Collins, S.R., Qu, H., Shales, M., Park, H.O., Hayles, J., et al., 2008. Conservation and rewiring of

functional modules revealed by an epistasis map in fission yeast. Science 322 (5900), 405–410. Sarasin-Filipowicz, M., Oakeley, E.J., Duong, F.H., Christen, V., Terracciano, L., Filipowicz, W., Heim, M.H., 2008. Interferon signaling and treatment outcome in chronic hepatitis C. Proc. Natl. Acad. Sci. U. S. A. 105 (19), 7034–7039. Shou, C., Bhardwaj, N., Lam, H.Y., Yan, K.K., Kim, P.M., Snyder, M., Gerstein, M.B., 2011. Measuring the evolutionary rewiring of biological networks. PLoS Comput. Biol. 7 (1), e1001050. Sim, S.C., Ingelman-Sundberg, M., 2011. Pharmacogenomic biomarkers: new tools in current and future drug therapy. Trends Pharmacol. Sci. 32 (2), 72–81. Sinha, R., Shaham, Y., Heilig, M., 2011. Translational and reverse translational research on the role of stress in drug craving and relapse. Psychopharmacology (Berl) 218 (1), 69–82. Subramanian, A., Tamayo, P., Mootha, V.K., Mukherjee, S., Ebert, B.L., Gillette, M.A., Paulovich, A., Pomeroy, S.L., Golub, T.R., Lander, E.S., et al., 2005. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. U. S. A. 102 (43), 15545–15550. Sun, M.G., Sikora, M., Costanzo, M., Boone, C., Kim, P.M., 2012. Network evolution: rewiring and signatures of conservation in signaling. PLoS Comput. Biol. 8 (3), e1002411. Tatebe, K., Zeytun, A., Ribeiro, R.M., Hoffmann, R., Harrod, K.S., Forst, C.V., 2010. Response network analysis of differential gene expression in human epithelial lung cells during avian influenza infections. BMC Bioinform. 11, 170. Taylor, M.W., Tsukahara, T., Brodsky, L., Schaley, J., Sanda, C., Stephens, M.J., McClintick, J.N., Edenberg, H.J., Li, L., Tavis, J.E., et al., 2007. Changes in gene expression during pegylated interferon and ribavirin therapy of chronic hepatitis C virus distinguish responders from nonresponders to antiviral therapy. J. Virol. 81 (7), 3391–3401. Taylor, M.W., Tsukahara, T., McClintick, J.N., Edenberg, H.J., Kwo, P., 2008. Cyclic changes in gene expression induced by Peg-interferon alfa-2b plus ribavirin in peripheral blood monocytes (PBMC) of hepatitis C patients during the first 10 weeks of treatment. J. Transl. Med. 6, 66. Torkamani, A., Schork, N.J., 2012. Background gene expression networks significantly enhance drug response prediction by transcriptional profiling. Pharmacogenomics J. 12 (5), 446–452. Vassilev, L.T., 2007. MDM2 inhibitors for cancer therapy. Trends Mol. Med. 13 (1), 23–31. Vaughan, L., 2001. Statistical Methods for the Information Professional. Information Today, Inc, Medford, NJ, USA. Wang, J., Huang, Q., Liu, Z.P., Wang, Y., Wu, L.Y., Chen, L., Zhang, X.S., 2011. NOA: a novel Network Ontology Analysis method. Nucleic Acids Res. 39 (13), e87. Wang, Y., Xu, M., Wang, Z., Tao, M., Zhu, J., Wang, L., Li, R., Berceli, S.A., Wu, R., 2012. How to cluster gene expression dynamics in response to environmental signals. Brief. Bioinform. 13 (2), 162–174. Welch, W.J., 1992. Mammalian stress response: cell physiology, structure/function of stress proteins, and implications for medicine and disease. Physiol. Rev. 72 (4), 1063–1081. Yu, X., Li, G., Chen, L., 2013. Prediction and early diagnosis of complex diseases by edge-network. Bioinformatics, http://dx.doi.org/10.1093/ bioinformatics/btt620. Yuasa, T., Takahashi, S., Hatake, K., Yonese, J., Fukui, I., 2011. Biomarkers to predict response to sunitinib therapy and prognosis in metastatic renal cell cancer. Cancer Sci. 102 (11), 1949–1957. Zeng, T., Chen, L., 2012. Tracing dynamic biological processes during phase transition. BMC Syst. Biol. 6 (Suppl. 1), S12. Zeng, T., Sun, S.Y., Wang, Y., Zhu, H., Chen, L., 2013. Network biomarkers reveal dysfunctional gene regulations during disease progression. FEBS J. 280 (22), 5682–5695. Zeng, T., Zhang, C.C., Zhang, W., Liu, R., Liu, J., Chen, L., 2014. Deciphering early development of complex diseases by progressive module network. Methods 67 (3), 334–343. Zhang, W., Zeng, T., Chen, L., 2014. EdgeMarker: identifying differentially correlated molecule pairs as edge-biomarkers. J. Theor. Biol., http://dx.doi.org/ 10.1016/j.jtbi.2014.05.041. Zheng, G., Tu, K., Yang, Q., Xiong, Y., Wei, C., Xie, L., Zhu, Y., Li, Y., 2008. ITFP: an integrated platform of mammalian transcription factors. Bioinformatics 24 (20), 2416–2417.