BACKGROUND
network-based approaches to the analysis of social and biological networks have provided scientists with powerful tools for the intuitive visualization of complex processes and for the principled integration of multiple data sources.
in genomics, network-based models – with genes represented by nodes and gene-gene interactions represented by edges – are inferred from experimental data or retrieved from manually curated repositories available online.
data-driven transcriptional networks, where a link between two nodes denotes a strong association between the corresponding gene expression profiles, are particularly useful as they represent a snapshot of the gene co-regulation in the experiment under study  <cit> .
within this context, scale free networks , in which the degree of connection of the member nodes follows a power law, have been widely used in the study of protein interactions  <cit> .
the construction of these networks usually relies on the computation of gene-gene correlations across replicate experiments, and on the subsequent thresholding of the absolute correlation values so as to define as connected only those genes with correlation above a chosen threshold  <cit> .
although this approach has proven extremely useful in identifying key hub genes in multiple biological conditions  <cit> , the high sensitivity of the obtained networks to the choice of threshold raises questions about the reproducibility of the obtained results, as well as about their biological meaning  <cit> .
an alternative approach is to use correlation networks , where all pairwise gene associations are considered, to avoid loss of information in those cases where the analysis focuses on the identification of groups of tightly connected genes , rather than on the identification of single key nodes .
while both types of graphs are referred to as networks in the following manuscript, it should be noted that the cn approach yields a fully interconnected  graph, on which some of the topological indices developed for “standard” networks with sparse links are not applicable.
regardless of the methodology used to infer the graph, network-derived gene modules can be investigated experimentally in order to gain insights into their biological function, or with the help of gene and pathway annotation resources.
additionally, the comparison of correlation networks from different conditions  may help identify modules whose connectivity is significantly altered in the compared conditions  <cit> .
connectivity-based comparisons may thus help identify “aggregate changes” that could be missed by standard methods of differential analysis comparing individual genes  <cit> .
in this study, we describe the development of a network-based analysis pipeline and its application to gene expression datasets from chemical perturbation experiments, with the goal of elucidating the modes of actions of the profiled perturbations.
we apply our pipeline to the analysis of the drugmatrix dataset from the national toxicology program   <cit>  and the tg-gates dataset from the japanese national institute of biomedical innovation  <cit> , two of the largest toxicogenomics datasets available, which contain organ-specific gene expression measurements for model organisms exposed to hundreds of chemical compounds with varying carcinogenicity and genotoxicity.
evidence accumulated to date has shown that machine learning techniques can successfully be applied to infer chemical carcinogenicity  from expression profiles of in vitro and in vivo assays  <cit> .
in our own previous work, we have shown that it is possible to infer highly accurate predictive models of chemical-associated long-term cancer risk from rat-based short-term toxicogenomics data, and to identify genes significantly associated with carcinogenesis  <cit> .
here, we aim to go beyond the inference of predictive models and the identification of single biomarker genes, towards the identification of gene modules or pathways significantly associated with the profiled chemical perturbations and the induced adverse phenotypes.
we do so by comparing the connectivity of gene modules in the networks derived from the control samples  to those obtained from samples collected after the exposure to specific chemical compounds.
to this end, we reconstruct chemical-specific transcriptional networks, and show that by grouping chemicals based on the similarity of their associated networks we can identify groups of chemicals or drugs with similar functions and similar carcinogenicity/genotoxicity profiles.
we also show that the in-silico annotation by pathway enrichment analysis of the gene modules with a differential connectivity  can point to the main molecular pathways induced by specific chemicals.
while network-based differential analysis of gene expression profiles is not new, the novelty of our approach lies in the application of the approach to aggregate chemical perturbations by their network similarity, in our application of a bootstrapping-based statistical significance test, and in the module-based meta-analysis of the differential connectivity results across multiple perturbations.
while in this study we focused on chemical-induced carcinogenesis and genotoxicity as the adverse phenotypes of interest, the approach we propose has broader applicability.
RESULTS
differential connectivity analysis of chemical perturbations
as shown in fig.
<dig>  the network-based pipeline for differential connectivity analysis starts by inferring chemical-specific compound networks, obtained from samples collected after the exposure to specific chemical compounds , and a network from the control samples, hereafter named “control network” .
groups of compounds are then identified based on the similarity of their individual network structures.
for each group, a new “aggregate compound network” is inferred by pooling all the samples across the clustered compounds .
next, modules of tightly connected genes are identified in each of the constructed networks and compared between conditions  in terms of module differential connectivity  .
given a module identified in one of the two networks under comparison , the mdc score is computed as the ratio of the average connectivity across all the genes within the module in the aggregate network  and in the control network .
this score represents changes of connectivity in the compound group with respect to the control .
mdc scores are computed for all modules in the networks, and tests of statistical significance and module specificity are performed to identify modules that will be further investigated through enrichment analysis based on pathway repositories and additional annotation sources .fig.
<dig> analysis workflow.
drugmatrix liver samples are used to infer chemical-specific compound networks  and a control network .
similarity in terms of network structure is evaluated to identify groups of compounds, whose samples are pooled to infer aggregate compound networks .
modules of tightly connected genes both in the control network and in compound aggregate networks are identified and compared across conditions in terms of module differential connectivity  .
modules with a change in connectivity that is highly specific to each compound group are investigated through pathway enrichment analysis
network comparison is sensitive to cross-dataset batch effects
before applying the described pipeline, we sought to rigorously assess the reproducibility of the methodologies used.
to this end, we carried out a systematic validation by comparing networks inferred from independent sets of samples profiling the same experimental conditions.
this approach was used under the assumption that networks derived from the same experimental conditions should in theory yield very similar structures, with differences in connectivity only due to sampling variability.
however, in practice many sources of “unwanted variation” can contribute to the measured differences, first among them differences due to non-reproducible experimental settings .
for this reason, we performed two types of validation experiments, aimed at identifying network inference and comparison methods showing a high reproducibility and robustness to possible batch effects.
first, networks derived from the control  liver samples of two independently generated datasets  were compared in terms of the composition of gene modules and their observed changes in connectivity.
for each of the modules identified from the tg-gates network, a correspondent module  was found in the drugmatrix network.
the high significance of the overlap for all the modules  confirmed that the two networks yielded the same structure and modules composition.
however, a systematic deviation of the distribution of differential connectivity values from one was detected in multiple modules.
a careful analysis showed that mdc values different from  <dig> were mainly caused by a difference in the distribution of correlation values between the drugmatrix and the tg-gates datasets .
furthermore, cross-dataset scaling or normalization could not eliminate these differences.
our evaluation points to the need for carefully assessing the methodologies and samples used for network inference and to adopt rigorous permutation-based statistical testing, in order to avoid interpreting artifacts or batch effects as actual differences in modules’ connectivity.
secondly, the inference and comparison methods were evaluated on networks derived from independent sample subsets extracted from the same dataset .
in this case, we attained more reproducible results, with the distribution of mdc values correctly centered at  <dig>  and with a lower variance when using the correlation network  approach rather than the scale-free network  approach .
an additional advantage of the simpler cn approach is that it does not require the selection of a threshold for the correlation values, a choice that is highly sensitive to the samples analyzed and strongly influences the subsequent calculation of mdc values.
based on these results, the subsequent analyses were all based on a single dataset , and the cn approach was selected as the network inference method of choice.
groups of similar chemicals can be inferred by network analysis
using the cn approach, we analyzed how different groups of compounds affect the connectivity pattern of specific gene modules.
first, we inferred the network from the non-treated liver samples, hereafter named “control network”.
the control network was clustered into  <dig> gene modules, with sizes ranging from  <dig> to  <dig> genes.
next, chemical-specific compound networks were inferred for each of the  <dig> chemicals for which at least ten replicate experiments  were available.
although a comparison of modules connectivity in control and individual chemicals-related networks was carried out , we aimed at identifying chemicals with similar network structure, which were then grouped to infer “aggregate compound networks”.
this step had two main goals: i) to study how well groups of chemicals with similar known features  can be identified through networks; and ii) to increase the sample size available for network inference.
for the aggregation of the chemicals, the compound networks were compared pairwise based on the similarity of their respective modules as measured by the adjusted rand index   <cit> .
the ari is a score specifically devised for the comparison of clustering results even when the two networks have different number of clusters .
hierarchical clustering was then applied to the matrix of ari’s to induce a similarity-based partial ordering and grouping of the chemicals .fig.
<dig> compounds aggregation.
similarity of  <dig> chemical compounds based on adjusted rand index .
a. heatmap of ari and grouping of compounds with similar networks structure.
b. zoom-in on the compounds grouping
the ari-based clustering yielded a clear separation between non-genotoxic carcinogens  and genotoxic non-carcinogens .
of notice, genotoxic compounds were not as well separated when we applied alternative, more standard clustering approaches, such as one based on the direct similarity of the chemicals’ expression profiles, and one based on the chemicals’ shared interacting proteins as annotated in public databases .
dynamic tree cutting  <cit>  was applied to the ari-based dendrogram with the parameters chosen so as to yield clusters with adequate sample size and sufficiently homogeneous sets of perturbations.
by this procedure,  <dig> chemical groups were identified and analyzed to assess their internal similarity as determined by multiple criteria.
according to drugbank annotation,  <dig> of the  <dig> chemical groups showed a high overlap in the main pharmacological action of their member chemicals as annotated in the stitch database  <cit>  .
as an example, group g <dig> includes only hypolipidemic drugs classified as “statins” present in the dataset, and it is separated from group g <dig>  which is mainly composed of “fibrates”, another class of antihyperlipidemic agents.
remarkably, some of these groups were not identified by standard gene expression-based analysis.
these include groups g <dig> and g <dig>  whose chemicals were separated into multiple clusters by dynamic tree cutting of the dendrogram shown in additional file 1: figure s <dig> a, obtained by gene expression similarity.
compounds within these two groups were indeed significantly similar, as measured by the overlap of their interacting proteins .
overall, a significant internal similarity was observed in  <dig> of the  <dig> groups for which protein annotations were available from public databases .
high overlap of shared side effects was also observed, although only two compound groups showed statistical significance by the permutation-based approach, due to a large number of side effects that are common to most of the compounds available in the sider database  <cit> .table  <dig> compound groups.
main functions retrieved through the stitch database
na = no annotations were found in ctd or in sider.
* = lower p-value cannot be obtained by chance
in summary, our analysis showed that the comparison of network modules by ari yields groups of chemicals with similar carcinogenicity-genotoxicity profiles in terms of pharmacological action, interacting proteins, and side effects.
module differential connectivity highlights chemicals’ modes of action
samples related to the  <dig> groups identified were used to infer aggregate compound networks, each representing the partial correlation among genes across all replicate experiments from the group of compounds considered.
as shown in fig.
<dig>  we aimed at comparing the control network with multiple perturbations by analyzing the changes in gene modules connectivity.
this analysis was repeated twice, first using the modules identified in the control network, and then using the modules identified in each of the aggregate compound networks.
control network-centered analysis
for each of the modules identified in the control network, the mdc score was computed to measure the change in connectivity  among the module’s genes due to the action of each group of chemicals.
significance of the mdc values was assessed by a bootstrap approach, whereby a confidence interval for each mdc is estimated by performing network inference on bootstrapped  versions of the original dataset, as detailed in methods and in additional file 1: figure s <dig>  an alternative, permutation-based significance testing procedure was also evaluated, whereby p-values are obtained by computing mdc values for random sets of genes, with sizes equal to each module under study.
however, this approach was not sufficiently stringent, yielding an excessively large number of likely false positive results .
the bootstrap approach identified a more parsimonious set of highly significant mdc values .
figure  <dig> summarizes the control network-centered results.
since the modules are defined in the control network, hence are the same for each pairwise comparison, the results can be represented as a matrix and an associated color-coded heatmap, with each row corresponding to a control network module, and each column corresponding to a compound group.
several modules manifest a remarkable change of connectivity, as captured by their mdc scores, and as confirmed by their estimated q-values .fig.
<dig> gain and loss of connectivity of control modules.
differential connectivity of  <dig> control modules  induced by  <dig> groups of compounds .
the heatmap is color-coded according to the mdc values, with blue and red indicating a loss and a gain of connectivity, respectively
of notice, the “turquoise” module was the only group showing a highly significant gain of connectivity for all the groups of compounds analyzed.
this result can be explained by the remarkably high number of genes included in this module , if compared to the median size of all the modules, equal to  <dig>  since chemical perturbations induced a global increase of the connectivity of the entire network , large size modules are more likely to show the same pattern, i.e.
a gain in the connectivity when compared to the control network.
moreover, biological annotation of this module showed a variety of different pathways, suggesting that compounds with diverse function could have an effect on this group of genes.
in order to focus on compound-specific effects, we computed a “specificity score” for each control module  <cit> .
the specificity score quantifies the uniqueness of a gain or loss of connectivity to a given compound group.
briefly, for a given group of compounds and a given module, the differences in mdc between that group of compounds and all the other groups are computed.
specificity is then defined as the sum of all the differences, with higher values identifying modules with a high mdc absolute value relative to all the others.
modules with scores exceeding a top percentile of the distribution of specificity values were subject to enrichment analysis and significant pathways were selected with fdr-corrected p-values.
additional file 3: table s <dig> shows enrichment analysis results of the best-ranked specific modules for each compound group, and a bipartite graph in fig.
<dig> is used to graphically represent the obtained associations between compound groups and enriched gene sets.
all groups except one  showed at least one top specific module significantly enriched for hallmarks gene sets .fig.
<dig> enrichment of specific control modules.
bipartite graph representing associations between compound groups and enriched hallmarks gene set corresponding to specifically altered modules extracted from the control network
aggregate network-centered analysis
taking an approach complementary to the one adopted in the control network-centered analysis, here a set of gene modules is defined for each of the aggregate compound networks.
that is, for each aggregate network, a set of densely connected modules is identified, and their change of connectivity with respect to the control network is calculated by mdc.
since a potentially distinct set of modules is identified in each aggregate network, this precludes the representation of the differential connectivity analysis results across the aggregate networks in matrix form.
in this analysis, we first identified high frequency  modules as those modules for which similar grouping of genes  was found across multiple aggregate networks .
a specificity score was then computed to highlight those with mdc values specific to a particular compound group.
we next identified, low frequency  modules, i.e., modules whose composition was unique to only one or few aggregate networks.
both hf and lf modules are reported in additional file 4: table s <dig>  and graphically represented in fig.
<dig>  where hallmarks gene sets have been used to investigate their biological function.
taken together, the findings described below confirm that the approach is capable of identifying known modes of action  <cit> , and of grouping compounds based on their coordinated effect on molecular pathways.fig.
<dig> enrichment of specific compounds-related modules.
bipartite graph representing associations between compound groups and enriched hallmark gene sets corresponding to specifically altered modules extracted from each aggregate compound network
in addition, a comparison of our network-based analysis results with those from standard differential expression analysis highlighted the complementarity of the two approaches.
differentially expressed genes identified as belonging to the “perturbational transcriptome” in our previous study  <cit> , and their correspondent enriched hallmarks, were evaluated in terms of their overlap with genes and gene sets identified by the present approach.
first, each of the control and compounds-related modules shown in additional file 3: table s <dig> and additional file 4: table s <dig> was scored for its enrichment in terms of differentially expressed genes.
as expected, a majority of the modules  were significantly enriched for differentially expressed genes .
additional file 5: table s <dig> shows a summary of this analysis where all genes contained in at least one of the modules identified are compared with the perturbational transcriptome.
despite the significant overlap yielded by the two approaches, a considerable number of genes were identified only by one of the methods, pointing to the complementarity of the approaches.
interestingly, many of the genes and hallmarks identified only with the network-based approach were associated to pathways previously implicated in mediation to chemical responses, and described in the following section, including heme metabolism, myc targets, and inflammation signaling.
discussion
the following sections aim to discuss the main biological findings resulting from both control-centered and compounds-centered analyses.
alcohol-induced liver inflammation
as shown in fig.
<dig>  g <dig> has a specific effect on a high number of pathways, mostly associated with inflammation and tumorigenesis.
given that liver samples have been considered in this study, the significant impact of this group on the entire pathway set could be explained by the high number of alcohols included in group g <dig>  in fact, alcohol-mediated activation of inflammation signaling pathways in the liver is known to increase tumorigenesis in mice and to activate pro-inflammatory cytokine, such as tumor necrosis factor alpha , interleukin  <dig> , and nuclear factor kappa b   <cit> .
this liver damage response has been reported for several members of the “g1-solvents” group, including allyl alcohol  <cit> , lipopolysaccaride, known as an endotoxin,  <cit> , and chloroform  <cit> .
in particular, chloroform has been associated to an anti-inflammatory action, possibly explaining the loss of connectivity of module “coral1” in additional file 4: table s <dig>  concordantly, tnf-α has been reported as showing a remarkable functional duality, being strongly engaged both in tissue regeneration/expansion and destruction  <cit> .
hypolipidemic compounds induce cholesterol metabolism and inflammation
groups g <dig>  and g <dig> , both including hypolipidemic compounds, show a specific alteration of modules related to cholesterol and fatty acid metabolism.
in particular, the highest specificity score in additional file 3: table s <dig> is obtained by the fibrates on module “thistle1”, enriched for fatty acid metabolism.
the mdc values color-coded in fig.
<dig> show that while all the other groups of chemicals cause a loc of the “thistle1” module, groups g <dig> and g <dig> are the only ones to produce a goc, with a higher mdc obtained by fibrates.
statins have a more significant effect on the module “honeydew”, enriched for cholesterol homeostasis , which has not a high specificity ranking for the fibrates, confirming different actions of those two distinct classes of drugs.
statins are also associated with stress response pathways, including oxidative phosphorylation, uv response and activation of tnf-α in both control  and compounds-related modules .
while hypoxia-inducible factors were found to play a role in the inhibition of cholesterol synthesis  <cit> , an inflammatory response of the liver has not been clearly reported in the literature.
however, the indication of liver damage as a rare side effect by fda might suggest that tissue-specific network analysis could capture most of the possible mechanisms induced by drugs exposure.
effect of estrogens, steroids and cancer drugs on cellular replication
both modules related to estrogens and to cancer treatment have an effect on the connectivity of module “maroon”, enriched for pathways related to cellular replication.
while this module is gaining connectivity as an effect of proliferation-inducing drugs contained in group g <dig> , a loss of connectivity is observed for g <dig>  and g <dig>  pointing to the disruption of the cellular replication machinery caused by anti-cancer drugs.
the gain of connectivity of g2m checkpoint induced by g <dig>  and p <dig> pathway in the compound-related modules  also confirms known mechanisms in the treatment of tumors and dna damage  <cit> .
the loss of connectivity of inflammation-related pathways observed on g <dig> perturbation could be explained by the known action of some steroids-related compounds.
in particular, a negative correlation with retinol and liver fibrosis has been previously reported  <cit> .
another interesting effect of one of the steroid compounds is revealed by enrichment analysis of modules altered by group g <dig>  mainly composed of antifungals.
our results show an effect of this group on coagulation pathways, observed in the control modules and supported by the literature  <cit> , as well as an effect on regeneration-related pathways, observed in the compounds-related modules.
remarkably, this group also contains mifepristone, a steroid compound that was also shown to have an effect on heme pathways.
in addition, significantly lower effects on coagulation pathways were shown for other steroid compounds with pharmacological action similar to mifepristone  <cit> .
this finding highlights a similar action of the compounds included in g <dig>  despite being assigned to different pharmacological classes.
non-homogeneous groups of compounds have known common effects
interesting results can be observed for non-homogeneous groups of compounds for which a predominant pharmacological action could not be assigned.
in particular, group g <dig>  shows a double effect, clearly visible in fig.
<dig>  by acting on inflammation-related pathways and on mitosis gene sets.
possible key players of this group are safrole, which has been shown to induce liver dna damage   <cit>  and methyl salicylate, shown to have an estrogenic potential  <cit> , thus increasing the coordinated activity of genes related to cellular replication pathways.
another interesting example is group g <dig>  containing four chemicals with apparently different functions.
as suggested by the specific loss of connectivity induced on module “coral1”, related to inflammatory response, the majority of compounds included in this group have antioxidant properties .
among these, atorvastatin is a compound belonging to the class of statins, which was not grouped with the other statins in g <dig>  interestingly, no other statins except for one were demonstrated to have antioxidant effects in vitro, confirming the grouping of compounds found with ari  <cit> .
methods
data sources
the drugmatrix  <cit> , available through the gene expression omnibus  with the accession number gse <dig>  contains gene expression profiles from male rat primary tissues  and cultured rat hepatocytes, corresponding to treatments with  <dig> chemicals, and including  <dig> control samples from rats kept in matched conditions.
each compound was administered at multiple doses and durations , and each combination of tissue, compound, time and dose was profiled in triplicate.
of the  <dig> chemicals tested,  <dig> were annotated with either carcinogenicity or genotoxicity information in the carcinogenic potency database   <cit> , corresponding to  <dig> profiles.
in our study, only the samples from liver were considered, both for controls  and for chemical perturbations represented by at least  <dig> samples and including all doses and durations available.
the average sample size used for network inference was of  <dig> samples, ranging from a minimum of  <dig> to a maximum of  <dig> samples, with the corresponding power to detect absolute correlation values of  <dig>  with significance below  <dig>  ranging from  <dig>  to  <dig> , respectively.
the toxicogenomics project-genomics assisted toxicity evaluation system   <cit> , is available through arrayexpress  and includes  <dig>  samples of male rat primary liver and kidney tissues, and cultured hepatocytes.
tg-gates tested  <dig> chemical compounds.
in our validation, only the samples from liver control were used, for a total of  <dig> samples used to infer a gene network.
data processing
both affymetrix datasets were normalized using the r bioconductor package frma and frmatools  <cit> .
the median absolute deviation  was used as the variation filter to select the  <dig> best-ranked probes whose expression was then considered for inference of transcriptional networks.
data normalization, gene selection, network inference and other analyses were performed with custom scripts developed using the programming languages r, and several bioconductor packages.
network inference and modules analysis
transcriptional network inference starts by defining an adjacency matrix a = {a
ij}, with weight a
ij denoting the strength of the relation of genes i and j in the expression data.
scale-free transformations  can then be applied to the correlation measurements to achieve a scale-free topology typical of biological networks, characterized by relatively few highly connected nodes  among a larger number of sparsely connected neighbors  <cit> .
in this work, we explored both the direct use of non-transformed correlation networks  as well as of scale-free transformed networks .
in order to obtain a correlation matrix, pearson correlation measures between all pairs of gene expression profiles were computed.
in the cn approach, the correlation matrix was directly used as the adjacency matrix a. conversely, in the sfn approach two additional steps were required: i) only those edges with correlation values exceeding a specific threshold were retained, with the threshold selected so that the resulting distribution of connectivities fitted a scale-free topology; and ii) the adjacency matrix a was computed by transforming the thresholded correlation values into a topological overlap matrix , which takes into account the indirect interactions between each couple of genes in the network  <cit> .
in both approaches, hierarchical clustering with ward’s method was then applied to the obtained adjacency matrix, and a dynamic tree cutting algorithm was applied to determine the number and composition of gene clusters, henceforth referred to as gene modules.
we used the r package cutreedynamic with minimum cluster size set to  <dig> genes, method set to “hybrid” and “deepsplit” parameter set to  <dig>  which allows a higher number of more homogeneous clusters if compared to other parameter settings.
finally, a module differential connectivity score  was used to compare the connectivity of gene modules between networks  <cit> .
for a specific module composed of n genes and an edge set en, the mdc measures the ratio of the weighted cardinalities of en in the two network, i.e.
:  <dig> mdcxy=enxeny
where |en
x| and |en
y| denote the average connectivities a
ijx and a
ijy among the module’s genes within networks x and y. mdc values below  <dig> represent a loss of connectivity of the module in x with respect to y, while values exceeding  <dig> indicate a gain of connectivity.
validation of inference methods
the reproducibility of the network inference procedure was evaluated by two approaches, under the assumption that networks derived from the same experimental conditions should yield very similar structures, with minimal differences in connectivity only due to sampling variability.
in the first validation, networks derived from the control  liver samples of two independent datasets, the drugmatrix  and the tg-gates , were constructed and compared.
the significance of the overlap between modules extracted from the two networks was assessed by means of the fisher exact test.
in the second validation, a resampling approach was adopted, whereby the control liver samples from the larger tg-gates dataset was randomly split into two equally sized datasets, and corresponding networks were derived and compared, with the procedure repeated  <dig> times.
in both validations, the distributions of mdcs obtained by the two inference approaches  were computed and compared.
inference of compound and aggregate compound networks
the inference approach  showing the best validation results was used to analyze how different compounds  affect the connectivity patterns of specific gene modules.
all non-treated liver samples available were used to construct the control network, while for the treatment-related networks we relied only on chemical compounds for which at least ten replicate experiments  were available, obtaining  <dig> compound networks.
in order to build aggregate networks representing multiple chemical compounds, the similarity of compound networks based on their module composition was measured by the adjusted rand index   <cit> .
the ari is a well-accepted measure that allows for the comparison of clustering results even when these yield different numbers of clusters.
groups of chemical compounds were then identified by applying dynamic tree cutting  to the hierarchical tree obtained by merging compounds, with the ward method, based on their ari similarity.
for each of the  <dig> groups detected, an aggregate network was built using partial correlation, in place of simple correlation, as the adjacency measure, so as to control for the potential confounding effect of the chemicals grouped.
the sample sizes of these groups ranged from  <dig> to  <dig>  with an average of  <dig> samples used to infer correlation values.
internal similarity of compounds in each group was assessed by evaluating the overlap of their interacting proteins, as retrieved through the ctd database  <cit> , by means of the fisher exact test.
specifically, a p-value was obtained for each pair of compounds in a group and the median p-value was used as the score of internal similarity among the aggregated set of compounds.
the significance of these measures was assessed by randomly selecting equally sized groups of compounds from the ctd database and computing their internal similarity, with the procedure repeated  <dig> times.
the same permutation-based approach was used to compute the significance of side effects similarity, this time considering sider  <cit>  as a knowledge source to retrieve compound-related information.
comparison with other clustering approaches
the grouping obtained with the network-based approach was compared to two alternative methodologies, based on: i) standard clustering methods which do not take advantage of the network structure; ii) information from available network-based knowledge sources.
the first approach starts by defining an average expression profile m
i
= , for each compound c
i, with values v
i obtained as mean expression values for each gene j across all samples related to c
i. pairwise distances between compounds c
i and c
k were thus estimated by computing the euclidean distance of the correspondent profiles m
i and m
k.
in the second approach, chemical-protein associations were retrieved through stitch, a widely adopted public repository of chemical-related networks  <cit> .
stitch contains multiple evidences of associations, with experimental evidences being one of the most reliable types of recorded interactions.
a network with all chemicals and their associated proteins was thus obtained by considering only associations coming from experimental evidences.
subsequently, a distance measure for compounds c
i and c
k was estimated as 1-j, j being the jaccard index of their set of associated proteins p
i and p
k,: j = p
i∩p
k
/p
i∪p
k
.
in both approaches, the pairwise distances computed were used as input to hierarchical clustering  to obtain compounds groupings.
selection and annotation of significant modules
for each control network-specific module, a confidence interval for the value of the corresponding mdc with respect to each aggregate network  was computed by randomly selecting with replacement the same number of samples from the replicates of non-treated samples.
after  <dig> iterations, the standard value of each log-transformed mdc was compared with the obtained estimates of mdc confidence interval.
the resulting p-values assess the significance of the deviation of the mdc from  <dig> .
the same bootstrap procedure was adopted for modules extracted from each aggregate network , this time by randomly selecting with replacement the same number of samples from the entire set of perturbations.
in order to rank the control modules most specifically altered by each compound group, the changes in connectivity of each module m measured for a compound group g
i with respect to the control network c were compared to those obtained for the other compound groups  <cit> .
first, the absolute value of the difference in the mdc scores between two groups of compounds gi and gj was computed for each module m as:  <dig> Δmgi,gj=logmdcmgi,c−logmdcmgj,c
this was used to compute the specificity of module m to compound group g
i, as:  <dig> spmgi=∑k=1nΔmgi,gk
where n denotes the total number of compound groups different from g
i. for each compound group, modules with score exceeding the top 5th percentile of the overall distribution of specificity scores were selected for enrichment analysis.
specificity of modules inferred from the aggregate compounds networks was assessed based on two alternative criteria, depending on the frequency of observation of the same module  in the networks.
modules identified in more than 50% of the aggregate networks were labeled as “high frequency” and a score spgi for module m to compound group g
i, was computed as described for the control network modules .
the specificity score was obtained as spmgintot where ntot denotes the total number of compound groups.
for each compound group, modules with score exceeding the top 5th percentile of the overall distribution of specificity scores were selected as “high frequency” specific modules.
modules identified in less than 50% of the networks were labelled as “low frequency” and selected based on significance of their correspondent mdc values .
both control network-derived and aggregate network-derived  specific modules were annotated by enrichment of the hallmark gene sets part of the msigdb compendium  <cit> .
significance of pathway enrichment was computed using a hyper-geometric distribution-based test and corrected for multiple hypothesis testing across multiple pathway gene sets via the false discovery rate  estimation.
hallmark pathways with fdr ≤  <dig>  was reported for each specific gene module in both control networks and aggregate networks.
selected genes and hallmark gene sets were compared with the “perturbational transcriptome”, a list of genes identified as significantly differentially expressed  in at least five compounds  <cit> .
each module was tested for enrichment of genes included in the perturbational transcriptome by a hyper-geometric test.
computational requirements
the pipeline was run on a linux node with two eight-core  <dig>  ghz intel xeon processors, and 128 gb ram.
several r scripts implementing the main steps illustrated in fig.
<dig> were run sequentially on the node.
network inference and modules identification  required approximately 17 s for the control network and 9 s for an individual compound network inferred from  <dig> samples.
pairwise similarity of  <dig> networks and grouping with the ari  took approximately 15 s. the run time for modules identification and differential connectivity analysis  between a 60-module control network and a 68-module aggregate compound network was about 6 s, with  <dig>  s needed for mdc computation of a module of  <dig> genes.
similar run times were required for each of the aggregate compound networks.
the most computationally intensive part of the pipeline was the significance testing based on bootstrapping, which required repeating correlation and mdc computation  <dig> times for each module .
considering an aggregate network inferred from  <dig> samples, this step took on average  <dig>  s/iteration, with a total of 28 min required for the complete assessment of mdc q-values, composed of  <dig> repetitions.
similarly, significance analysis of compounds-related modules took  <dig>  s/iteration, with a total run time of 30 min needed for q-values estimation.
in summary, the control network-centered analysis took on average  <dig>  min and the aggregate network-centered analysis for a single group of chemicals took an average of  <dig>  min.
it should be noted that once the control and aggregate networks are inferred, the statistical significance tests for the control network-centered and each of the aggregate network-centered analyses are independent of each other and can thus be executed in parallel.
CONCLUSIONS
we have presented a pipeline for transcriptional network inference and comparison that was primarily designed for the analysis of chemical perturbations from high-throughput transcriptional screening experiments.
here, we applied it to the analysis of gene expression profiles from rat-based chemical exposure experiments.
we show that groups of chemicals with similar functions and carcinogenicity/genotoxicity profiles can be identified through our proposed pipeline.
in addition, modules with altered connectivity due to the action of specific compounds were enriched for pathways actually related to the chemicals’ action.
these findings highlight potential advantages in the application of this network-based approach.
in the context of drug discovery , the methods presented here could help assign new functions to novel  drugs, based on the similarity of their associated network with those built for other known compounds.
additionally, networks with patients as nodes could be compared with the same tools in order to identify groups with a similar response to a set of drugs.
in fact, the proposed methodology has broad applicability beyond the uses here described and could be used as an alternative or as a complement to standard approaches of differential gene expression analysis.
as currently implemented, the pipeline identifies modules with differential connectivity irrespective of link directionality.
an extension of the method possibly worth investigating would be the modeling of directional links, e.g., obtained by distinguishing between positively and negatively correlated gene pairs.
this would require a modification of the connectivity score here adopted, and could potentially highlight modules for which a specific perturbation is able to change the direction of nodes connectivity with respect to a reference condition.
on the other hand, the obtained results would likely become more difficult to interpret and annotate.
additional files
additional file 1:
supplementary figures.
additional figures with supplementary results.
additional file 2: table s <dig>  global differential connectivity  of the entire networks inferred for compound groups with respect to the control network.
additional file 3: table s <dig>  control network-related modules with connectivity specifically altered by compound groups.
additional file 4: table s <dig>  aggregate network-related modules with connectivity specifically altered by compound groups.
additional file 5: table s <dig>  comparison of genes and gene sets identified as differentially connected in the network-based approach and by the standard differential expression analysis.
abbreviations
ariadjusted rand index
cncorrelation network
mdcmodule differential connectivity
sfnscale free network

 network-based approaches to the analysis of social and biological networks have provided scientists with powerful tools for the intuitive visualization of complex processes and for the principled integration of multiple data sources. in genomics, network-based models – with genes represented by nodes and gene-gene interactions represented by edges – are inferred from experimental data or retrieved from manually curated repositories available online. data-driven transcriptional networks, where a link between two nodes denotes a strong association between the corresponding gene expression profiles, are particularly useful as they represent a snapshot of the gene co-regulation in the experiment under study. within this context, scale free networks, in which the degree of connection of the member nodes follows a power law, have been widely used in the study of protein interactions. the construction of these networks usually relies on the computation of gene-gene correlations across replicate experiments, and on the subsequent thresholding of the absolute correlation values so as to define as connected only those genes with correlation above a chosen threshold. although this approach has proven extremely useful in identifying key hub genes in multiple biological conditions, the high sensitivity of the obtained networks to the choice of threshold raises questions about the reproducibility of the obtained results, as well as about their biological meaning. an alternative approach is to use correlation networks, where all pairwise gene associations are considered, to avoid loss of information in those cases where the analysis focuses on the identification of groups of tightly connected genes, rather than on the identification of single key nodes. while both types of graphs are referred to as networks in the following manuscript, it should be noted that the cn approach yields a fully interconnected graph, on which some of the topological indices developed for “standard” networks with sparse links are not applicable. regardless of the methodology used to infer the graph, network-derived gene modules can be investigated experimentally in order to gain insights into their biological function, or with the help of gene and pathway annotation resources. additionally, the comparison of correlation networks from different conditions may help identify modules whose connectivity is significantly altered in the compared conditions. connectivity-based comparisons may thus help identify “aggregate changes” that could be missed by standard methods of differential analysis comparing individual genes. in this study, we describe the development of a network-based analysis pipeline and its application to gene expression datasets from chemical perturbation experiments, with the goal of elucidating the modes of actions of the profiled perturbations. we apply our pipeline to the analysis of the drugmatrix dataset from the national toxicology program and the tg-gates dataset from the japanese national institute of biomedical innovation, two of the largest toxicogenomics datasets available, which contain organ-specific gene expression measurements for model organisms exposed to hundreds of chemical compounds with varying carcinogenicity and genotoxicity. evidence accumulated to date has shown that machine learning techniques can successfully be applied to infer chemical carcinogenicity from expression profiles of in vitro and in vivo assays. in our own previous work, we have shown that it is possible to infer highly accurate predictive models of chemical-associated long-term cancer risk from rat-based short-term toxicogenomics data, and to identify genes significantly associated with carcinogenesis. here, we aim to go beyond the inference of predictive models and the identification of single biomarker genes, towards the identification of gene modules or pathways significantly associated with the profiled chemical perturbations and the induced adverse phenotypes. we do so by comparing the connectivity of gene modules in the networks derived from the control samples to those obtained from samples collected after the exposure to specific chemical compounds. to this end, we reconstruct chemical-specific transcriptional networks, and show that by grouping chemicals based on the similarity of their associated networks we can identify groups of chemicals or drugs with similar functions and similar carcinogenicity/genotoxicity profiles. we also show that the in-silico annotation by pathway enrichment analysis of the gene modules with a differential connectivity can point to the main molecular pathways induced by specific chemicals. while network-based differential analysis of gene expression profiles is not new, the novelty of our approach lies in the application of the approach to aggregate chemical perturbations by their network similarity, in our application of a bootstrapping-based statistical significance test, and in the module-based meta-analysis of the differential connectivity results across multiple perturbations. while in this study we focused on chemical-induced carcinogenesis and genotoxicity as the adverse phenotypes of interest, the approach we propose has broader applicability.  differential connectivity analysis of chemical perturbations as shown in fig.  the network-based pipeline for differential connectivity analysis starts by inferring chemical-specific compound networks, obtained from samples collected after the exposure to specific chemical compounds, and a network from the control samples, hereafter named “control network”. groups of compounds are then identified based on the similarity of their individual network structures. for each group, a new “aggregate compound network” is inferred by pooling all the samples across the clustered compounds. next, modules of tightly connected genes are identified in each of the constructed networks and compared between conditions in terms of module differential connectivity. given a module identified in one of the two networks under comparison, the mdc score is computed as the ratio of the average connectivity across all the genes within the module in the aggregate network and in the control network. this score represents changes of connectivity in the compound group with respect to the control. mdc scores are computed for all modules in the networks, and tests of statistical significance and module specificity are performed to identify modules that will be further investigated through enrichment analysis based on pathway repositories and additional annotation sources.fig. analysis workflow. drugmatrix liver samples are used to infer chemical-specific compound networks and a control network. similarity in terms of network structure is evaluated to identify groups of compounds, whose samples are pooled to infer aggregate compound networks. modules of tightly connected genes both in the control network and in compound aggregate networks are identified and compared across conditions in terms of module differential connectivity. modules with a change in connectivity that is highly specific to each compound group are investigated through pathway enrichment analysis network comparison is sensitive to cross-dataset batch effects before applying the described pipeline, we sought to rigorously assess the reproducibility of the methodologies used. to this end, we carried out a systematic validation by comparing networks inferred from independent sets of samples profiling the same experimental conditions. this approach was used under the assumption that networks derived from the same experimental conditions should in theory yield very similar structures, with differences in connectivity only due to sampling variability. however, in practice many sources of “unwanted variation” can contribute to the measured differences, first among them differences due to non-reproducible experimental settings. for this reason, we performed two types of validation experiments, aimed at identifying network inference and comparison methods showing a high reproducibility and robustness to possible batch effects. first, networks derived from the control liver samples of two independently generated datasets were compared in terms of the composition of gene modules and their observed changes in connectivity. for each of the modules identified from the tg-gates network, a correspondent module was found in the drugmatrix network. the high significance of the overlap for all the modules confirmed that the two networks yielded the same structure and modules composition. however, a systematic deviation of the distribution of differential connectivity values from one was detected in multiple modules. a careful analysis showed that mdc values different from were mainly caused by a difference in the distribution of correlation values between the drugmatrix and the tg-gates datasets. furthermore, cross-dataset scaling or normalization could not eliminate these differences. our evaluation points to the need for carefully assessing the methodologies and samples used for network inference and to adopt rigorous permutation-based statistical testing, in order to avoid interpreting artifacts or batch effects as actual differences in modules’ connectivity. secondly, the inference and comparison methods were evaluated on networks derived from independent sample subsets extracted from the same dataset. in this case, we attained more reproducible results, with the distribution of mdc values correctly centered at and with a lower variance when using the correlation network approach rather than the scale-free network approach. an additional advantage of the simpler cn approach is that it does not require the selection of a threshold for the correlation values, a choice that is highly sensitive to the samples analyzed and strongly influences the subsequent calculation of mdc values. based on these results, the subsequent analyses were all based on a single dataset, and the cn approach was selected as the network inference method of choice. groups of similar chemicals can be inferred by network analysis using the cn approach, we analyzed how different groups of compounds affect the connectivity pattern of specific gene modules. first, we inferred the network from the non-treated liver samples, hereafter named “control network”. the control network was clustered into gene modules, with sizes ranging from to genes. next, chemical-specific compound networks were inferred for each of the chemicals for which at least ten replicate experiments were available. although a comparison of modules connectivity in control and individual chemicals-related networks was carried out, we aimed at identifying chemicals with similar network structure, which were then grouped to infer “aggregate compound networks”. this step had two main goals: i) to study how well groups of chemicals with similar known features can be identified through networks; and ii) to increase the sample size available for network inference. for the aggregation of the chemicals, the compound networks were compared pairwise based on the similarity of their respective modules as measured by the adjusted rand index. the ari is a score specifically devised for the comparison of clustering results even when the two networks have different number of clusters. hierarchical clustering was then applied to the matrix of ari’s to induce a similarity-based partial ordering and grouping of the chemicals.fig. compounds aggregation. similarity of chemical compounds based on adjusted rand index. a. heatmap of ari and grouping of compounds with similar networks structure. b. zoom-in on the compounds grouping the ari-based clustering yielded a clear separation between non-genotoxic carcinogens and genotoxic non-carcinogens. of notice, genotoxic compounds were not as well separated when we applied alternative, more standard clustering approaches, such as one based on the direct similarity of the chemicals’ expression profiles, and one based on the chemicals’ shared interacting proteins as annotated in public databases. dynamic tree cutting was applied to the ari-based dendrogram with the parameters chosen so as to yield clusters with adequate sample size and sufficiently homogeneous sets of perturbations. by this procedure, chemical groups were identified and analyzed to assess their internal similarity as determined by multiple criteria. according to drugbank annotation, of the chemical groups showed a high overlap in the main pharmacological action of their member chemicals as annotated in the stitch database. as an example, group g includes only hypolipidemic drugs classified as “statins” present in the dataset, and it is separated from group g which is mainly composed of “fibrates”, another class of antihyperlipidemic agents. remarkably, some of these groups were not identified by standard gene expression-based analysis. these include groups g and g whose chemicals were separated into multiple clusters by dynamic tree cutting of the dendrogram shown in additional file 1: figure s a, obtained by gene expression similarity. compounds within these two groups were indeed significantly similar, as measured by the overlap of their interacting proteins. overall, a significant internal similarity was observed in of the groups for which protein annotations were available from public databases. high overlap of shared side effects was also observed, although only two compound groups showed statistical significance by the permutation-based approach, due to a large number of side effects that are common to most of the compounds available in the sider database.table compound groups. main functions retrieved through the stitch database na = no annotations were found in ctd or in sider. * = lower p-value cannot be obtained by chance in summary, our analysis showed that the comparison of network modules by ari yields groups of chemicals with similar carcinogenicity-genotoxicity profiles in terms of pharmacological action, interacting proteins, and side effects. module differential connectivity highlights chemicals’ modes of action samples related to the groups identified were used to infer aggregate compound networks, each representing the partial correlation among genes across all replicate experiments from the group of compounds considered. as shown in fig.  we aimed at comparing the control network with multiple perturbations by analyzing the changes in gene modules connectivity. this analysis was repeated twice, first using the modules identified in the control network, and then using the modules identified in each of the aggregate compound networks. control network-centered analysis for each of the modules identified in the control network, the mdc score was computed to measure the change in connectivity among the module’s genes due to the action of each group of chemicals. significance of the mdc values was assessed by a bootstrap approach, whereby a confidence interval for each mdc is estimated by performing network inference on bootstrapped versions of the original dataset, as detailed in methods and in additional file 1: figure s an alternative, permutation-based significance testing procedure was also evaluated, whereby p-values are obtained by computing mdc values for random sets of genes, with sizes equal to each module under study. however, this approach was not sufficiently stringent, yielding an excessively large number of likely false positive results. the bootstrap approach identified a more parsimonious set of highly significant mdc values. figure  summarizes the control network-centered results. since the modules are defined in the control network, hence are the same for each pairwise comparison, the results can be represented as a matrix and an associated color-coded heatmap, with each row corresponding to a control network module, and each column corresponding to a compound group. several modules manifest a remarkable change of connectivity, as captured by their mdc scores, and as confirmed by their estimated q-values.fig. gain and loss of connectivity of control modules. differential connectivity of control modules induced by groups of compounds. the heatmap is color-coded according to the mdc values, with blue and red indicating a loss and a gain of connectivity, respectively of notice, the “turquoise” module was the only group showing a highly significant gain of connectivity for all the groups of compounds analyzed. this result can be explained by the remarkably high number of genes included in this module, if compared to the median size of all the modules, equal to since chemical perturbations induced a global increase of the connectivity of the entire network, large size modules are more likely to show the same pattern, i.e. a gain in the connectivity when compared to the control network. moreover, biological annotation of this module showed a variety of different pathways, suggesting that compounds with diverse function could have an effect on this group of genes. in order to focus on compound-specific effects, we computed a “specificity score” for each control module. the specificity score quantifies the uniqueness of a gain or loss of connectivity to a given compound group. briefly, for a given group of compounds and a given module, the differences in mdc between that group of compounds and all the other groups are computed. specificity is then defined as the sum of all the differences, with higher values identifying modules with a high mdc absolute value relative to all the others. modules with scores exceeding a top percentile of the distribution of specificity values were subject to enrichment analysis and significant pathways were selected with fdr-corrected p-values. additional file 3: table s shows enrichment analysis results of the best-ranked specific modules for each compound group, and a bipartite graph in fig.  is used to graphically represent the obtained associations between compound groups and enriched gene sets. all groups except one showed at least one top specific module significantly enriched for hallmarks gene sets.fig. enrichment of specific control modules. bipartite graph representing associations between compound groups and enriched hallmarks gene set corresponding to specifically altered modules extracted from the control network aggregate network-centered analysis taking an approach complementary to the one adopted in the control network-centered analysis, here a set of gene modules is defined for each of the aggregate compound networks. that is, for each aggregate network, a set of densely connected modules is identified, and their change of connectivity with respect to the control network is calculated by mdc. since a potentially distinct set of modules is identified in each aggregate network, this precludes the representation of the differential connectivity analysis results across the aggregate networks in matrix form. in this analysis, we first identified high frequency modules as those modules for which similar grouping of genes was found across multiple aggregate networks. a specificity score was then computed to highlight those with mdc values specific to a particular compound group. we next identified, low frequency modules, i.e., modules whose composition was unique to only one or few aggregate networks. both hf and lf modules are reported in additional file 4: table s and graphically represented in fig.  where hallmarks gene sets have been used to investigate their biological function. taken together, the findings described below confirm that the approach is capable of identifying known modes of action, and of grouping compounds based on their coordinated effect on molecular pathways.fig. enrichment of specific compounds-related modules. bipartite graph representing associations between compound groups and enriched hallmark gene sets corresponding to specifically altered modules extracted from each aggregate compound network in addition, a comparison of our network-based analysis results with those from standard differential expression analysis highlighted the complementarity of the two approaches. differentially expressed genes identified as belonging to the “perturbational transcriptome” in our previous study, and their correspondent enriched hallmarks, were evaluated in terms of their overlap with genes and gene sets identified by the present approach. first, each of the control and compounds-related modules shown in additional file 3: table s and additional file 4: table s was scored for its enrichment in terms of differentially expressed genes. as expected, a majority of the modules were significantly enriched for differentially expressed genes. additional file 5: table s shows a summary of this analysis where all genes contained in at least one of the modules identified are compared with the perturbational transcriptome. despite the significant overlap yielded by the two approaches, a considerable number of genes were identified only by one of the methods, pointing to the complementarity of the approaches. interestingly, many of the genes and hallmarks identified only with the network-based approach were associated to pathways previously implicated in mediation to chemical responses, and described in the following section, including heme metabolism, myc targets, and inflammation signaling. discussion the following sections aim to discuss the main biological findings resulting from both control-centered and compounds-centered analyses. alcohol-induced liver inflammation as shown in fig.  g has a specific effect on a high number of pathways, mostly associated with inflammation and tumorigenesis. given that liver samples have been considered in this study, the significant impact of this group on the entire pathway set could be explained by the high number of alcohols included in group g in fact, alcohol-mediated activation of inflammation signaling pathways in the liver is known to increase tumorigenesis in mice and to activate pro-inflammatory cytokine, such as tumor necrosis factor alpha, interleukin, and nuclear factor kappa b. this liver damage response has been reported for several members of the “g1-solvents” group, including allyl alcohol, lipopolysaccaride, known as an endotoxin,, and chloroform. in particular, chloroform has been associated to an anti-inflammatory action, possibly explaining the loss of connectivity of module “coral1” in additional file 4: table s concordantly, tnf-α has been reported as showing a remarkable functional duality, being strongly engaged both in tissue regeneration/expansion and destruction. hypolipidemic compounds induce cholesterol metabolism and inflammation groups g and g, both including hypolipidemic compounds, show a specific alteration of modules related to cholesterol and fatty acid metabolism. in particular, the highest specificity score in additional file 3: table s is obtained by the fibrates on module “thistle1”, enriched for fatty acid metabolism. the mdc values color-coded in fig.  show that while all the other groups of chemicals cause a loc of the “thistle1” module, groups g and g are the only ones to produce a goc, with a higher mdc obtained by fibrates. statins have a more significant effect on the module “honeydew”, enriched for cholesterol homeostasis, which has not a high specificity ranking for the fibrates, confirming different actions of those two distinct classes of drugs. statins are also associated with stress response pathways, including oxidative phosphorylation, uv response and activation of tnf-α in both control and compounds-related modules. while hypoxia-inducible factors were found to play a role in the inhibition of cholesterol synthesis, an inflammatory response of the liver has not been clearly reported in the literature. however, the indication of liver damage as a rare side effect by fda might suggest that tissue-specific network analysis could capture most of the possible mechanisms induced by drugs exposure. effect of estrogens, steroids and cancer drugs on cellular replication both modules related to estrogens and to cancer treatment have an effect on the connectivity of module “maroon”, enriched for pathways related to cellular replication. while this module is gaining connectivity as an effect of proliferation-inducing drugs contained in group g, a loss of connectivity is observed for g and g pointing to the disruption of the cellular replication machinery caused by anti-cancer drugs. the gain of connectivity of g2m checkpoint induced by g and p pathway in the compound-related modules also confirms known mechanisms in the treatment of tumors and dna damage. the loss of connectivity of inflammation-related pathways observed on g perturbation could be explained by the known action of some steroids-related compounds. in particular, a negative correlation with retinol and liver fibrosis has been previously reported. another interesting effect of one of the steroid compounds is revealed by enrichment analysis of modules altered by group g mainly composed of antifungals. our results show an effect of this group on coagulation pathways, observed in the control modules and supported by the literature, as well as an effect on regeneration-related pathways, observed in the compounds-related modules. remarkably, this group also contains mifepristone, a steroid compound that was also shown to have an effect on heme pathways. in addition, significantly lower effects on coagulation pathways were shown for other steroid compounds with pharmacological action similar to mifepristone. this finding highlights a similar action of the compounds included in g despite being assigned to different pharmacological classes. non-homogeneous groups of compounds have known common effects interesting results can be observed for non-homogeneous groups of compounds for which a predominant pharmacological action could not be assigned. in particular, group g shows a double effect, clearly visible in fig.  by acting on inflammation-related pathways and on mitosis gene sets. possible key players of this group are safrole, which has been shown to induce liver dna damage and methyl salicylate, shown to have an estrogenic potential, thus increasing the coordinated activity of genes related to cellular replication pathways. another interesting example is group g containing four chemicals with apparently different functions. as suggested by the specific loss of connectivity induced on module “coral1”, related to inflammatory response, the majority of compounds included in this group have antioxidant properties. among these, atorvastatin is a compound belonging to the class of statins, which was not grouped with the other statins in g interestingly, no other statins except for one were demonstrated to have antioxidant effects in vitro, confirming the grouping of compounds found with ari. methods data sources the drugmatrix, available through the gene expression omnibus with the accession number gse contains gene expression profiles from male rat primary tissues and cultured rat hepatocytes, corresponding to treatments with chemicals, and including control samples from rats kept in matched conditions. each compound was administered at multiple doses and durations, and each combination of tissue, compound, time and dose was profiled in triplicate. of the chemicals tested, were annotated with either carcinogenicity or genotoxicity information in the carcinogenic potency database, corresponding to profiles. in our study, only the samples from liver were considered, both for controls and for chemical perturbations represented by at least samples and including all doses and durations available. the average sample size used for network inference was of samples, ranging from a minimum of to a maximum of samples, with the corresponding power to detect absolute correlation values of with significance below ranging from to, respectively. the toxicogenomics project-genomics assisted toxicity evaluation system, is available through arrayexpress and includes samples of male rat primary liver and kidney tissues, and cultured hepatocytes. tg-gates tested chemical compounds. in our validation, only the samples from liver control were used, for a total of samples used to infer a gene network. data processing both affymetrix datasets were normalized using the r bioconductor package frma and frmatools. the median absolute deviation was used as the variation filter to select the best-ranked probes whose expression was then considered for inference of transcriptional networks. data normalization, gene selection, network inference and other analyses were performed with custom scripts developed using the programming languages r, and several bioconductor packages. network inference and modules analysis transcriptional network inference starts by defining an adjacency matrix a = {a ij}, with weight a ij denoting the strength of the relation of genes i and j in the expression data. scale-free transformations can then be applied to the correlation measurements to achieve a scale-free topology typical of biological networks, characterized by relatively few highly connected nodes among a larger number of sparsely connected neighbors. in this work, we explored both the direct use of non-transformed correlation networks as well as of scale-free transformed networks. in order to obtain a correlation matrix, pearson correlation measures between all pairs of gene expression profiles were computed. in the cn approach, the correlation matrix was directly used as the adjacency matrix a. conversely, in the sfn approach two additional steps were required: i) only those edges with correlation values exceeding a specific threshold were retained, with the threshold selected so that the resulting distribution of connectivities fitted a scale-free topology; and ii) the adjacency matrix a was computed by transforming the thresholded correlation values into a topological overlap matrix, which takes into account the indirect interactions between each couple of genes in the network. in both approaches, hierarchical clustering with ward’s method was then applied to the obtained adjacency matrix, and a dynamic tree cutting algorithm was applied to determine the number and composition of gene clusters, henceforth referred to as gene modules. we used the r package cutreedynamic with minimum cluster size set to genes, method set to “hybrid” and “deepsplit” parameter set to which allows a higher number of more homogeneous clusters if compared to other parameter settings. finally, a module differential connectivity score was used to compare the connectivity of gene modules between networks. for a specific module composed of n genes and an edge set en, the mdc measures the ratio of the weighted cardinalities of en in the two network, i.e.: mdcxy=enxeny where |en x| and |en y| denote the average connectivities a ijx and a ijy among the module’s genes within networks x and y. mdc values below represent a loss of connectivity of the module in x with respect to y, while values exceeding indicate a gain of connectivity. validation of inference methods the reproducibility of the network inference procedure was evaluated by two approaches, under the assumption that networks derived from the same experimental conditions should yield very similar structures, with minimal differences in connectivity only due to sampling variability. in the first validation, networks derived from the control liver samples of two independent datasets, the drugmatrix and the tg-gates, were constructed and compared. the significance of the overlap between modules extracted from the two networks was assessed by means of the fisher exact test. in the second validation, a resampling approach was adopted, whereby the control liver samples from the larger tg-gates dataset was randomly split into two equally sized datasets, and corresponding networks were derived and compared, with the procedure repeated times. in both validations, the distributions of mdcs obtained by the two inference approaches were computed and compared. inference of compound and aggregate compound networks the inference approach showing the best validation results was used to analyze how different compounds affect the connectivity patterns of specific gene modules. all non-treated liver samples available were used to construct the control network, while for the treatment-related networks we relied only on chemical compounds for which at least ten replicate experiments were available, obtaining compound networks. in order to build aggregate networks representing multiple chemical compounds, the similarity of compound networks based on their module composition was measured by the adjusted rand index. the ari is a well-accepted measure that allows for the comparison of clustering results even when these yield different numbers of clusters. groups of chemical compounds were then identified by applying dynamic tree cutting to the hierarchical tree obtained by merging compounds, with the ward method, based on their ari similarity. for each of the groups detected, an aggregate network was built using partial correlation, in place of simple correlation, as the adjacency measure, so as to control for the potential confounding effect of the chemicals grouped. the sample sizes of these groups ranged from to with an average of samples used to infer correlation values. internal similarity of compounds in each group was assessed by evaluating the overlap of their interacting proteins, as retrieved through the ctd database, by means of the fisher exact test. specifically, a p-value was obtained for each pair of compounds in a group and the median p-value was used as the score of internal similarity among the aggregated set of compounds. the significance of these measures was assessed by randomly selecting equally sized groups of compounds from the ctd database and computing their internal similarity, with the procedure repeated times. the same permutation-based approach was used to compute the significance of side effects similarity, this time considering sider as a knowledge source to retrieve compound-related information. comparison with other clustering approaches the grouping obtained with the network-based approach was compared to two alternative methodologies, based on: i) standard clustering methods which do not take advantage of the network structure; ii) information from available network-based knowledge sources. the first approach starts by defining an average expression profile m i  =, for each compound c i, with values v i obtained as mean expression values for each gene j across all samples related to c i. pairwise distances between compounds c i and c k were thus estimated by computing the euclidean distance of the correspondent profiles m i and m k. in the second approach, chemical-protein associations were retrieved through stitch, a widely adopted public repository of chemical-related networks. stitch contains multiple evidences of associations, with experimental evidences being one of the most reliable types of recorded interactions. a network with all chemicals and their associated proteins was thus obtained by considering only associations coming from experimental evidences. subsequently, a distance measure for compounds c i and c k was estimated as 1-j, j being the jaccard index of their set of associated proteins p i and p k,: j = p i∩p k /p i∪p k. in both approaches, the pairwise distances computed were used as input to hierarchical clustering to obtain compounds groupings. selection and annotation of significant modules for each control network-specific module, a confidence interval for the value of the corresponding mdc with respect to each aggregate network was computed by randomly selecting with replacement the same number of samples from the replicates of non-treated samples. after iterations, the standard value of each log-transformed mdc was compared with the obtained estimates of mdc confidence interval. the resulting p-values assess the significance of the deviation of the mdc from. the same bootstrap procedure was adopted for modules extracted from each aggregate network, this time by randomly selecting with replacement the same number of samples from the entire set of perturbations. in order to rank the control modules most specifically altered by each compound group, the changes in connectivity of each module m measured for a compound group g i with respect to the control network c were compared to those obtained for the other compound groups. first, the absolute value of the difference in the mdc scores between two groups of compounds gi and gj was computed for each module m as: Δmgi,gj=logmdcmgi,c−logmdcmgj,c this was used to compute the specificity of module m to compound group g i, as: spmgi=∑k=1nΔmgi,gk where n denotes the total number of compound groups different from g i. for each compound group, modules with score exceeding the top 5th percentile of the overall distribution of specificity scores were selected for enrichment analysis. specificity of modules inferred from the aggregate compounds networks was assessed based on two alternative criteria, depending on the frequency of observation of the same module in the networks. modules identified in more than 50% of the aggregate networks were labeled as “high frequency” and a score spgi for module m to compound group g i, was computed as described for the control network modules. the specificity score was obtained as spmgintot where ntot denotes the total number of compound groups. for each compound group, modules with score exceeding the top 5th percentile of the overall distribution of specificity scores were selected as “high frequency” specific modules. modules identified in less than 50% of the networks were labelled as “low frequency” and selected based on significance of their correspondent mdc values. both control network-derived and aggregate network-derived specific modules were annotated by enrichment of the hallmark gene sets part of the msigdb compendium. significance of pathway enrichment was computed using a hyper-geometric distribution-based test and corrected for multiple hypothesis testing across multiple pathway gene sets via the false discovery rate estimation. hallmark pathways with fdr ≤  was reported for each specific gene module in both control networks and aggregate networks. selected genes and hallmark gene sets were compared with the “perturbational transcriptome”, a list of genes identified as significantly differentially expressed in at least five compounds. each module was tested for enrichment of genes included in the perturbational transcriptome by a hyper-geometric test. computational requirements the pipeline was run on a linux node with two eight-core  ghz intel xeon processors, and 128 gb ram. several r scripts implementing the main steps illustrated in fig.  were run sequentially on the node. network inference and modules identification required approximately 17 s for the control network and 9 s for an individual compound network inferred from samples. pairwise similarity of networks and grouping with the ari took approximately 15 s. the run time for modules identification and differential connectivity analysis between a 60-module control network and a 68-module aggregate compound network was about 6 s, with  s needed for mdc computation of a module of genes. similar run times were required for each of the aggregate compound networks. the most computationally intensive part of the pipeline was the significance testing based on bootstrapping, which required repeating correlation and mdc computation times for each module. considering an aggregate network inferred from samples, this step took on average  s/iteration, with a total of 28 min required for the complete assessment of mdc q-values, composed of repetitions. similarly, significance analysis of compounds-related modules took  s/iteration, with a total run time of 30 min needed for q-values estimation. in summary, the control network-centered analysis took on average  min and the aggregate network-centered analysis for a single group of chemicals took an average of  min. it should be noted that once the control and aggregate networks are inferred, the statistical significance tests for the control network-centered and each of the aggregate network-centered analyses are independent of each other and can thus be executed in parallel.  we have presented a pipeline for transcriptional network inference and comparison that was primarily designed for the analysis of chemical perturbations from high-throughput transcriptional screening experiments. here, we applied it to the analysis of gene expression profiles from rat-based chemical exposure experiments. we show that groups of chemicals with similar functions and carcinogenicity/genotoxicity profiles can be identified through our proposed pipeline. in addition, modules with altered connectivity due to the action of specific compounds were enriched for pathways actually related to the chemicals’ action. these findings highlight potential advantages in the application of this network-based approach. in the context of drug discovery, the methods presented here could help assign new functions to novel drugs, based on the similarity of their associated network with those built for other known compounds. additionally, networks with patients as nodes could be compared with the same tools in order to identify groups with a similar response to a set of drugs. in fact, the proposed methodology has broad applicability beyond the uses here described and could be used as an alternative or as a complement to standard approaches of differential gene expression analysis. as currently implemented, the pipeline identifies modules with differential connectivity irrespective of link directionality. an extension of the method possibly worth investigating would be the modeling of directional links, e.g., obtained by distinguishing between positively and negatively correlated gene pairs. this would require a modification of the connectivity score here adopted, and could potentially highlight modules for which a specific perturbation is able to change the direction of nodes connectivity with respect to a reference condition. on the other hand, the obtained results would likely become more difficult to interpret and annotate. additional files additional file 1: supplementary figures. additional figures with supplementary results. additional file 2: table s global differential connectivity of the entire networks inferred for compound groups with respect to the control network. additional file 3: table s control network-related modules with connectivity specifically altered by compound groups. additional file 4: table s aggregate network-related modules with connectivity specifically altered by compound groups. additional file 5: table s comparison of genes and gene sets identified as differentially connected in the network-based approach and by the standard differential expression analysis. abbreviations ariadjusted rand index cncorrelation network mdcmodule differential connectivity sfnscale free network 
additionally, the method has broad applicability beyond the uses here described and could be used as an alternative or as a complement to standard approaches of differential gene expression analysis. starting from chemical-specific transcriptional networks inferred from these data, we show that the proposed comparative analysis of their associated networks identifies groups of chemicals with similar functions and similar carcinogenicity/genotoxicity profiles. methods for inference and comparison of biological networks are emerging as powerful tools for the identification of groups of tightly connected genes whose activity may be altered during disease progression or due to chemical perturbations. we also show that the in-silico annotation by pathway enrichment analysis of the gene modules with a significant gain or loss of connectivity for specific groups of compounds can reveal molecular pathways significantly associated with the chemical perturbations and their likely modes of action. correlation networkstoxicogenomicscompounds similaritygene expressioncomparative analysischemical perturbationsjoint 14th network tools and applications in biology workshop and 11th integrative bioinformatics symposium nettab/ib 2015bari, italy 14- october issue-copyright-statement© the author 2017 the proposed pipeline for transcriptional network inference and comparison is highly reproducible and allows grouping chemicals with similar functions and carcinogenicity/genotoxicity profiles. in this study, we describe a pipeline for network comparison and its application to the analysis of gene expression datasets from chemical perturbation experiments, with the goal of elucidating the modes of actions of the profiled perturbations. we apply our pipeline to the analysis of the drugmatrix and the tg-gates, two of the largest toxicogenomics resources available, containing gene expression measurements for model organisms exposed to hundreds of chemical compounds with varying carcinogenicity and genotoxicity. in the context of drug discovery or drug repositioning, the methods presented here could help assign new functions to novel or existing drugs, based on the similarity of their associated network with those built for other known compounds.
 methods for inference and comparison of biological networks are emerging as powerful tools for the identification of groups of tightly connected genes whose activity may be altered during disease progression or due to chemical perturbations. connectivity-based comparisons help identify aggregate changes that would be difficult to detect with differential analysis methods comparing individual genes. methods in this study, we describe a pipeline for network comparison and its application to the analysis of gene expression datasets from chemical perturbation experiments, with the goal of elucidating the modes of actions of the profiled perturbations. we apply our pipeline to the analysis of the drugmatrix and the tg-gates, two of the largest toxicogenomics resources available, containing gene expression measurements for model organisms exposed to hundreds of chemical compounds with varying carcinogenicity and genotoxicity.  starting from chemical-specific transcriptional networks inferred from these data, we show that the proposed comparative analysis of their associated networks identifies groups of chemicals with similar functions and similar carcinogenicity/genotoxicity profiles. we also show that the in-silico annotation by pathway enrichment analysis of the gene modules with a significant gain or loss of connectivity for specific groups of compounds can reveal molecular pathways significantly associated with the chemical perturbations and their likely modes of action.  the proposed pipeline for transcriptional network inference and comparison is highly reproducible and allows grouping chemicals with similar functions and carcinogenicity/genotoxicity profiles. in the context of drug discovery or drug repositioning, the methods presented here could help assign new functions to novel or existing drugs, based on the similarity of their associated network with those built for other known compounds. additionally, the method has broad applicability beyond the uses here described and could be used as an alternative or as a complement to standard approaches of differential gene expression analysis. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords correlation networkstoxicogenomicscompounds similaritygene expressioncomparative analysischemical perturbationsjoint 14th network tools and applications in biology workshop and 11th integrative bioinformatics symposium nettab/ib 2015bari, italy 14- october issue-copyright-statement© the author 2017 