Pathway analysis

Once the differentially expressed genes have been determined, the most common analytical method is pathway analysis to figure out which ‘biological functions’ are altered after compound exposure (Fig. 2d and e). However, the definition of ‘pathways’ or ‘biological functions’ depends on the database used. These definitions are evolving and may even be considered as somewhat arbitrary.89 Pathways are species dependent and so care must be taken when using the databases to ensure that the appropriate organism-specific pathway or ontology databases are available. The two most common analytical methods to determine the differentially expressed pathways or functions are the simple hypergeometric enrichment test, and Gene Set Enrichment Analysis (GSEA, Fig. 2d and e).90 The difference between the two lies in the null hypothesis. The hypergeometric enrichment test (reviewed recently91) investigates whether a pathway or a biological function occurs more often in DEGs compared to an appropriate background: usually either the set of genes measured in the microarray/RNA-Seq experiment or the entire genome of the species in the database (Fig. 2d). The null hypothesis is that the genes of a pathway are not enriched in the DEGs. Therefore, this method requires a predefined cut-off for determining which genes are significantly differentially expressed (see previously). GSEA, on the other hand, uses the expression value of all measured genes. It ranks the genes according to a metric (e.g. fold change) and then determines whether the genes from a set (e.g. from a pathway) occur in the high or low end of the ranked list. The null hypothesis here is that the genes from the set occur randomly in the ranked list. GSEA uses a Kolmogorov–Smirnov test for statistical significance of the enrichment. GSEA does not require a pre-defined cut-off to be specified for DEGs, in contrast to simple enrichment analysis (Fig. 2e).90

Both methods require gene sets for testing. Such gene sets can be obtained from the different pathway databases available, some of which are summarized in Table 3. Many of the toxicogenomic studies mentioned earlier have used pathway analysis, illustrating how this can be carried out using different pathway databases after selection of DEGs. Gene Ontology74,99 (GO) is probably the most commonly used database.100–102 However its hierarchical nature, as well as the nonspecificity of the higher GO layers, lead to difficulty in interpretation. To avoid such difficulties, it is good practice according to the authors' experience to group the annotations and identify the common grounds of the found GO terms; this feature is found in many online GO enrichment tools.

All pathway-based enrichment methods have a curation bias: the most important genes or pathways are well researched so they tend to have more ontological entries, or in the case of pathways, more member genes. Enrichment analysis by default does not give entirely novel mechanisms of action because some understanding about the genes involved needs to be provided to annotate them with meaningful pathways. However, the method contextualises experimental findings with the currently available biological insight. It is a common problem to receive a large number of ‘enriched’ pathways from such analysis, so the choice of background correction and filtering for relevant mechanisms is frequently employed. Kim et al. used the GOrilla103 tool to examine altered pathways after MPP+ induction in human neuroblastoma cells, finding different nucleosome assembly Gene Ontology biological processes to be enriched in a time-dependent manner.74

After enriched pathways have been found, a pathway map can be formed to shed light on causative biological events. An example is the work of Bell et al.,104 where the authors conducted a DEG analysis by determining FC from TG-GATEs.18 These DEGs were used to calculate the enriched pathways in Reactome constructing a “computationally predicted adverse outcome pathway” for each compound for a specific pathological phenotype. The usefulness of the method was validated with the example of the fatty liver disease caused by carbon tetrachloride.

Abdul Hameed et al. used protein interaction networks (see relevant section below) and pathway analysis to find toxic pathways involved in liver injuries in rats in vivo.101 They showed decreased metabolism in the liver but increased inflammatory pathway activity and increased expression of genes in fibrosis-relevant pathways (Fig. 3 – replica from ref. 101).

Open in a separate window

Fig. 3

An example network of enriched KEGG pathways of compound-induced differentially expressed genes relating to liver fibrosis. Of the 15 pathways shown, down-regulated pathways were predominantly metabolism related; whereas up-regulated pathways were related to processes associated with liver fibrosis, for example, the focal adhesion pathway is pathway 3 (P3) and also immune related pathways are depicted (P7-antigen processing or P9-chemokine signaling). The authors hypothesised that the metabolic pathways could be related to external factors (e.g. altered food intake) or an indication of reduced liver function. Genes with an average fold change >1.5 are in red, <0.75 are in green, and the remaining in grey. PX represents pathway number X and is represented in a blue circle. DOI: ; 10.1371/journal.pone.0112193.g003.

Pathway analysis forms the basis of most toxicogenomics analyses. The results of it may not be trivial to understand; determining mode-of-action from hundreds of DEGs and hundreds of enriched pathways is not always possible. As such, further methods have been developed to annotate experimental gene expression data with additional context.

Compound signature matching methods

A further development in toxicogenomic methods is signature-matching approaches, where compound-induced gene expression signatures are evaluated against a pre-existing compound signature library in order to make predictions about their potential toxicity. Compound signature matching methods have been used in a broad range of applications from side effect prediction19 to drug repurposing,105 based on the assumption that compounds inducing similar gene expression signatures will have similar effects in a biological system. In the field of toxicity, this allows the matching of test compounds to those with known toxicity profiles, or that have a known mechanism of toxicity. Importantly, basing the comparison on transcriptomic read-outs, rather than compound structure, may lead to a similarity profile very distinct from that obtained by structural similarity.87

Transcriptomic profiles of compounds can be obtained from compound signature collections such as CMap19 or LINCS (Table 1).51,52 As these collections measure compound-induced gene expression in vitro, a greater number of compounds can be queried when compared to the in vivo measurements in the toxicity-specific databases mentioned above such as TG-GATEs or DrugMatrix. In this part of the review, we therefore will focus on the in vitro databases CMap and LINCS, and their utilization in understanding and predicting compound toxicity.

Using these signature libraries, researchers can measure the similarity between compounds in gene expression space. A widely-used method to do this is connectivity mapping,19 which takes into account that the most strongly differentially expressed genes are likely to be more informative than the entire transcriptome. Connectivity mapping describes the enrichment of a ‘query’ signature (for instance, a list of the top most up- and down-regulated genes) against a reference transcriptomic profile (e.g. of a known toxicant) (Fig. 4). This is measured by a connectivity score based on the Kolmogorov–Smirnov statistic for the up- and down-regulated genes of a query compound. The original paper describing CMap illustrated how connectivity mapping could be used to elucidate the mechanism of action of a compound or predict side effects such as weight gain,19 and several early applications of connectivity mapping in toxicology are covered in a mini-review by Smalley et al.106

Open in a separate window

Fig. 4

Connectivity mapping for compound signatures. Lists of up- and down-regulated genes resulting from perturbation by a compound are compared against gene expression signatures from reference compounds. Positive connectivity (where the genes up-regulated by the compound under test are also up-regulated by the reference compound, and similarly for down-regulated genes) indicates that the two compounds induce similar gene expression profiles; negative connectivity indicates the opposite.

More recently, a case study of the use of gene expression data in drug discovery projects described how such an approach was used to evaluate the toxicity of compounds inhibiting PDE10A, an antipsychotic target.107 Expression profiling was carried out on human embryonic kidney (HEK293) cells for the compounds under development, revealing a strong downregulation of tubulin genes. The level of tubulin downregulation correlated with high levels of micronuclei formation, suggesting that the tubulin genes could be used as a predictive signature of micronuclei formation. These signature genes were then used to query the Connectivity Map to find compounds with similar patterns of gene expression. Four of the five most similar compounds returned by this approach were known genotoxic compounds, one of which is commonly used as a positive reference in the micronuclei formation test. This result was used to suggest subsequent transcriptomic experiments, which validated the link between the tubulin genes and micronuclei formation. The authors suggest that transcriptomic profiling could therefore provide an early indicator of potential genotoxicity, allowing compounds to be excluded well before the micronucleus test, which is usually performed late in the drug development pipeline.107

As well as testing for connectivity to known toxic compounds, compound signature similarity can be used to infer mechanisms of toxicity.108 One case study involves the use of connectivity to predict novel hERG (human ether-a-go-go-related gene) K+ channel inhibitors.88 Inhibition of the hERG channel leads to an increased risk of sudden cardiac death,109 but known hERG inhibitors are diverse with respect to their structure and primary targets, causing difficulty in the computational identification of potential inhibitory compounds.88 In order to investigate whether transcriptomic signatures could provide a signal of hERG inhibition, CMap profiles of 673 drugs including 119 known hERG inhibitors were clustered using affinity propagation, a clustering algorithm based on the idea of communication between data points.110 Similarities in the profiles of structurally diverse known hERG inhibitors were used to create a transcriptomic profile of hERG inhibition in different cell lines, revealing differential expression in groups of genes enriched for diverse processes including cholesterol and isoprenoid biosynthesis and the cell cycle. Clusters enriched for hERG inhibitors predicted novel inhibitors that showed significantly greater inhibition than randomly selected compounds, illustrating how CMap data can be used to generate signatures of toxicity based entirely on public data.

As well as the general issues faced in the analysis of gene expression data (as described above), there are further considerations arising from the use of in vitro cell line measurements in the largest compound-induced signature databases, CMap and LINCS. It is known that gene expression in cell lines does not always correlate closely with that measured in the corresponding organ;100 further, the gene expression response to compounds may be affected by the type of cell line used.108 Differences in cell line response, as well as between dosage and time point of compound administration, must therefore be taken into account when analysing this type of compound-induced signature. Nonetheless, as demonstrated in this section, signature-matching approaches can be a powerful tool for early hypothesis generation before later in vivo validation.

Utilizing biological networks for toxicogenomic studies

Biological interaction networks, such as protein–protein interaction or signaling networks, can be useful tools to decipher the mechanism of toxicity. Biological networks can be directed when we know which way the information flows from one node (protein, miRNA, small molecule, gene, etc.) to the other; or undirected when this information is unknown or it has no meaning, e.g. proteins forming a complex. Biological networks, especially directed signaling networks, allow us to follow the cellular response of a compound treatment from the compound's target to the differentially expressed genes. Different biological networks and databases are compiled from various data sources with varying coverage and information content, e.g. whether a network is directed or whether an interaction is inhibitory or excitatory (signed) (Fig. 5(i) and Table 4). The most commonly used biological networks are protein–protein interaction (PPI) networks, whose nodes represent proteins and edges represent interactions i.e. the binding of one protein to another. The researcher in every toxicogenomics project has to determine whether they want to look into a specific toxicological process deeply or map a general response and choose the network accordingly.

Open in a separate window

Fig. 5

Inferring mode of toxicity using network biology. The initial step (i) is to select a network with appropriate coverage and information content according to the question at hand. Subsequently, (ii) gene expression data needs to be merged with the selected PPI network database. Following this, (iii) an algorithm connects the differentially expressed genes in the network. The resulting putative toxicity networks can be depicted and suggest the mode of toxicity for a compound (iv in yellow), but further experimental validation is required to confirm the prediction. DEG – differentially expressed gene, ACO – ant colony optimization. For the different databases, see Table 4.

Table 4

Network resources for toxicogenomic studies

Network resource name	Description	Number of interactions in human	Species	Web address
Biomodels111	Small-scale dataset containing rate-related interactions. Varying coverage by model type.	Varies in scale 10–1000	Various	https://www.ebi.ac.uk/biomodels-main
NRF2Ome114	Small scale manually curated oxidative stress and NRF2 response specific database. Interactions are directed and signed.	289 NRF2 specific PPI	Human	http://nrf2.elte.hu
OmniPath89	Manually curated partly directed and signed signaling database which integrates other high quality interaction sources.	50247	Human	http://omnipathdb.or/
SignaLink2115	Multilayer signaling database with regulations and predicted interactions. The manually curated interactions are directed and signed.	1640 high confidence PPI	Human, fruit fly, C. elegans	signalink.org
Signor116	Manually curated pathway interactions with directions and signs.	19312	Human, mouse, rat	signor.uniroma2.it
Reactome117	Manually curated large scale reaction centered pathway database, which focuses on protein complexes, but the interactions are directed and signed.	11426	Different organisms including, human, rat, mouse	www.reactome.org
HPRD118	Historic, no longer updated database of manually curated undirected interactions.	41327	Human	www.hprd.org
Bioplex119	Large-scale immunopurification and mass spectrometry based protein interaction database.	70000	Human	http://bioplex.hms.harvard.edu
BioGRID120,121	Genetic and protein interactions from low and high throughput publications.	406487	Multiple species including human, rat, mouse	https://thebiogrid.org
IntAct120,121	Large scale protein interaction database collection.	310183	Mostly human, but contains other species data as well including mouse	https://www.ebi.ac.uk/intact
InWeb_IM122	Large scale collection of PPI datasets with orthological predictions.	625640	Human	https://www.intomics.com/inbio/map/#home
HAPPI-2123	Large database collection of protein interactions with a confidence score.	2922202	Human	http://discovery.informatics.uab.edu/HAPPI
mentha124	Scored collection of interactions from publications and databases.	309088	Multiple model organisms including mouse, rat, and humans	http://mentha.uniroma2.it
STRING125	Large scale predicted and curated interactions database. Uses text mining and orthology to cover the interactions in different species.	11353056	Multiple different organisms	https://string-db.org
Yu et al.126	Inferred high-confidence human protein–protein interactions from multiple data sources.	80980	Human interactions	https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-13-79

Open in a separate window

The interactions from biological networks can be characterized based on the source of the interaction and the types of annotation available, such as the direction, strength, kinetics, and sign (inhibitory or activating) of an interaction. While the ultimate aim for biological network studies is to model the whole cell and organism using detailed quantitative interactions, as of yet such detailed models are only available for a few genes or proteins in the BioModels database,111–113 rendering them unsuitable for toxicogenomics modelling at the current stage.

Manually curated databases typically contain somewhat less detailed information. Such databases include HPRD118 for undirected human interactions, and OmniPath89 or Signor116 for directed and signed signaling information. Reactome117 assembles pathways from curated interactions, but in some cases the directionality is impossible to define. Such manually curated databases are biased toward well-studied proteins and interactions, but these tend to be more accurate than high throughput databases. On the other hand, some databases contain information obtained from large, high-throughput experiments, such as BioGRID,120 BIOPLEX,119 and MINT.121 Although these databases contain many interactions, not every interaction is manually checked, so the confidence is usually lower. Most of the experiments are derived from yeast two hybrid model systems, which do not cover nuclear interactions or interactions in the cell membrane; an exception is BIOPLEX, which uses immunoaffinity purification with mass-spectrometry which gives unbiased, reliable data, but cannot differentiate the exact formation of complexes. However, the advantage of such high throughput databases is that they provide unbiased and large-scale information.

Other interaction databases aim to aggregate information from multiple sources, such as STRING, InWeb_IM and HAPPI.122,125,127 Appropriate filtering of such databases can make them applicable to answer toxicological questions, but it should be noted that merging different databases can increase the noise as well as the coverage. Consistency of data and annotations from multiple sources is a frequently recurring problem in this case.

To utilise biological networks (which are chosen according to relevant criteria, as in Fig. 5 step i), DEGs or transcriptomic signatures are first matched to proteins (Fig. 5 step ii). Identifier matching tools, like the UniProt retrieve128 service or the Protein Identifier Cross-Reference resource,129 can help to do this step.

The next step is to identify which functions these proteins affect in the network (Fig. 5 step iii). Most methods use random walk with restart algorithms, including ENRICHNET,130 NETPEA,131 and NetWalk.132 A related approach is the heat diffusion based algorithms such as HotNet2133 or DMFIND.134 Random walk with restart begins from the protein equivalents of the selected DEGs and walks around the PPI graph, with a random chance of restarting, to see which proteins can be reached from the start. In a case study testing the NetWalk algorithm, Komurov et al. used a unified PPI and transcription factor-target gene network.132 They captured the cell cycle arresting function of p53 to sublethal doses of doxorubicin and the apoptosis induction of p53 at lethal doses in MCF7 cell lines. HotNet2 was developed and successfully used for module assignment in pan-cancer data. It detected 16 such modules including the p53 and the NOTCH signaling module in multiple cancers.133

A more sophisticated method to find the affected proteins in the network is the Ant Colony Optimization (ACO),135 where the random walker (ant) leaves a ‘pheromone trail’ behind it, which increases the probability that the next ant will walk the same path. The strength of the pheromone trail depends on a function of the visited nodes in the graph. For example, if an ant reaches another signature node – such as a DEG – then the next ant can walk the same path and connect the new signature nodes with a path. It is an extension of random walk methods because ACO can connect, in the network sense, distant paths and not just discover the neighbourhood of the signature nodes.

In the toxicogenomics field, ACO was successfully used by Abdul Hameed et al.101 to uncover how toxicants can cause liver fibrosis through extracellular matrix bound growth factors. The authors determined differentially expressed genes using the rank product method from DrugMatrix17 data and also clustered them based on their co-expression in liver fibrosis. The differentially expressed genes and the co-expressed genes from the enriched clusters were mapped to a previously inferred and rescored high quality PPI network.126 KeyPathwayMiner,136 an ACO implementation for network analysis, was next used to construct the liver fibrosis-associated network. The network was then clustered with the EAGLE137 algorithm implemented in the Clusterviz Cytoscape plugin138 to find the network module most highly correlated with liver fibrosis. This method was shown to uncover novel interactions in liver fibrosis, which could not be revealed using pathway enrichment of co-expressed and differentially expressed genes. In this module, the extracellular matrix compartments and bounded growth factors were overrepresented, which was validated via independent data sets. The utilization of a PPI network enhanced the scope of the analysis, because it incorporated the indirectly affected genes, whose expression themselves was unchanged.

With network biology tools, the feedback effects can be followed from the targets of toxicants to the measured gene expression signature through transcription factors and a putative adverse outcome pathway can be constructed. Melas et al.102 achieved that in the case of drug-induced lung injury. They used Reactome as a source of protein interactions and a collection of transcription factor target data to connect gene expression signatures of drugs from CMAP with transcription factors. This analysis used a modified Integer Linear Programming algorithm.139 Integer Linear Programming is a tree-growing algorithm that finds the shortest tree between two sets of nodes in a directed graph. In this study, the algorithm was modified by adding transcription factors as a third set of nodes that have to be reached. The trees start to grow from the targets of toxicants, through transcription factors, to the differentially expressed genes. These trees formed the putative adverse outcome pathways for specific compounds. They validated their method with an independent pathway growing algorithm and random controls. The developed trees identified central apoptosis relevant proteins such as p53, CASP3, BCL2, BAX, CASP6, CASP8, CASP9 etc. and key signaling proteins such as FOS and JUN. Furthermore, these paths showed potential targets to avoid drug-induced lung injury and the authors tested specified drugs, which counteracted the lung injury as a toxic endpoint.

Biological networks can help to uncover hidden modes of toxicity with the help of gene expression data. They work with the assumption that the level of a transcript's expression is highly correlated with the amount of protein. However, this assumption is not absolutely true in all cases.140 To choose a proper biological network for a toxicogenomics study, the researcher must choose between information content and coverage. Nonspecific interaction databases with large coverage are suitable to generate unbiased hypotheses in toxicogenomics. If the coverage is not so important but the information content and reliability is a key issue, then manually curated database such as Signor or OmniPath or even smaller databases like NRF2Ome114 may be more appropriate. If kinetic modelling is the aim then the researcher must initially look up a relevant model from the BioModels database. The middle ground could be a database such as Reactome to show specific toxic responses with an appropriate coverage of signaling in humans and model organisms.

The authors acknowledge the useful comments from Fredrik Svensson and others from Bender Group. We thank Eszter Ari's help suggesting papers. BAD would like to acknowledge Tim James for useful discussions. LLP acknowledges Lorentz Jäntschi for useful suggestions and support. DM is funded by European Research Council grant number 336159. BAD is funded by the EPSRC (grant number 1827220). EO is funded by the BBSRC (grant number 1501561).