Introduction

Therapeutic drugs, industrial and environmental chemicals, and other compounds can have undesired toxic effects in diverse parts of the human body. This might be due to their direct interaction with cell components, but also, if they are metabolized, due to effects of their downstream catabolism [1]. These direct and indirect effects depend on the tissues considered (e.g., hepatotoxicity in liver [2] or nephrotoxicity in kidney [3]), on the human individual genetic variation (e.g., if a catabolizing enzyme has variable specificity associated to a single nucleotide polymorphism [4]) and on environmental factors (e.g., interaction with nutrients [5]). Thus, the toxic effects a given compound exerts in a given organ are a result of very complex interactions at cellular, organismic and environmental level, involving often multiple toxicity mechanisms [6].

The organ-specific study of drug effects on human physiology, while necessary, is however difficult given the impossibility of systematically testing the human system. Therefore, the commonly established approaches in the field of drug development and toxicology rely on studies of drug effects in animal organisms (in vivo) or in cultivated animal and human cell lines (in vitro). Profiling of gene expression can inform at a very detailed level about the effects of drugs. With the advent of high throughput technologies (in particular, the cost-effective microarrays), genome-wide transcription profiles of specific drug effects have been captured across multiple organs or cell types, drug dose levels and drug treatment time intervals. This enabled toxicologists to study the biochemical and biological mechanisms affected by a specific drug as well as the mechanisms behind organ-specific and general drug-induced toxicity. Most importantly, large-scale initiatives (e.g., DrugMatrix by Iconix Biosciences [7], the Japanese toxicogenomics project [8] and the EU funded InnoMed PredTox project [9]) have produced highly controlled and standardized toxicological experiments with rodent organisms (rat and mouse) and rodent and human cell lines, for hundreds of chemicals, resulting in comprehensive toxicogenomics resources for systematic analysis of drug induced toxicity effects, such as DrugMatrix [7], CMap [10] and Open TG-GATEs [11].

In this respect, the aim of the recently established field of toxicogenomics is to integrate these rapidly expanding collections of data (including gene expression profiles, clinical chemistry measurements, histopathology assessment data and curated literature annotations) in order to identify more specific and sensitive human toxicity biomarkers that should lead to more accurate toxicity detection at early stage of drug development and ultimately to safe marketed drugs [12].

While their contribution to our understanding of drug-induced toxicity is undeniable, one fundamental question that is difficult to address comprehensibly regarding the different systems used for experimentation in toxicogenomics studies is how reliably the specific animal organism, model system (in vivo or in vitro), reflects the toxicity effects of drugs for a given human organ in terms of genes (or molecular functions) whose expression is affected. In particular, it would be desirable to have a protocol to evaluate the appropriateness of multiple model systems so that the optimal one can be chosen for each organ, depending on the ability of the given model to affect the expression of genes and functions of interest.

Here, we describe a method to obtain lists of ranked genes according to gene expression changes in a toxicogenomics dataset profiling gene expression changes upon multiple drugs, conditions, and time points. If the dataset uses a particular tissue or cellular line as model to test drug toxicity, the expectation is that the ranked list reflects genes generally associated to toxicity in this model. Because our goal is to model toxicity in human organs, if the model is non-human, we take the corresponding human orthologous genes. Then, we optimize the methodology and parameters to obtain this ranking by comparison to the set of genes already known to be associated to drug toxicity in the corresponding organ in the literature.

Due to the multi-factorial experimental designs of the existing toxicological databases, we need to integrate data across multiple model systems, drugs, dose levels and time points. A relatively simple way is to summarize the genome-wide drug-induced transcriptional changes in prioritized (ranked) gene lists, followed by a rank-aggregation approach [13] to combine them into a single organ-specific gene list. While this approach is commonly used for prioritization of disease-related genes [14], here we use it as baseline method for prioritization of toxicity-related genes. In particular, we exploit the robustness (with respect to noise, outliers and errors in the data to be aggregated) and computational efficiency of recently proposed probabilistic methods for rank-aggregation, i.e., RRA [15] and BIRRA [16].

We applied this approach to a toxicogenomics dataset corresponding to 33 compounds (mainly toxic drugs) profiled in primary human hepatocytes, primary rat hepatocytes, rat liver (as models for hepatotoxicity) and rat kidney (as model for nephrotoxicity), obtained from the Open TG-GATEs database. The latter was especially suitable for our purposes due to the consistent experimental design applied across the four model systems. We evaluated multiple ways to aggregate these data for the description of drug-induced toxicity in human kidney and liver.

In doing this analysis, we discovered some genes that while ranking very high and very often in response to liver and kidney drug toxicity, have not been yet mentioned in the literature in relation to toxicity. These could be potentially working in novel toxicity-related pathways. We discuss the striking differences observed between human and rat hepatocytes, and the insight that can be gained into which processes can be modelled when using cell lines (rat hepatocytes) versus entire organs (rat liver). We use these results to describe the suitability of each of these models to study different molecular functions associated to organ toxicity according to the sets of genes whose expression is affected.

Methods

Unless otherwise stated, the complete protocol for data analysis was performed using the R version 3.2 computing environment [17] and packages from the open-source software development project Bioconductor version 3.1 [18].

Results

To develop a protocol for gene prioritization according to organ-specific toxicity from large-scale panels of toxicogenomics profiles, we integrated microarray gene expression time series from the Open TG-GATEs database [11] corresponding to 33 compounds (the majority of which are nephrotoxic or hepatotoxic drugs, S1 File) for four model systems: human hepatocytes (HH), rat hepatocytes (RH), rat liver (RL) and rat kidney (RK). We chose to analyze HH, RH, and RL to be able to compare multiple models for hepatotoxicity (in vivo versus in vitro, and in human versus rat), and then added RK as model for nephrotoxicity, to be able to compare results for two different in vivo models (liver versus kidney).

After microarray preprocessing and gene expression modeling, we performed a differential analysis to identify genes with significant differential expression profiles over time for each compound and model system. This step was crucial for the protocol as it provides the compound-specific gene ranks that serve as input for the rank-aggregation step.

Using different rank-aggregation methods, compound-wise gene rankings were summarized into one ranked gene list, which was limited to one-to-one orthologs between human and rat (see Methods, subsection Gene ranking for details). We evaluated the performance of the rank-aggregation methods against human liver and kidney toxicity-related genes derived from co-occurrence analysis of PubMed abstract annotations (see Methods, subsection Toxicity-related genes from the literature for details). We further compared the rank-aggregation methods based on the receiver-operating characteristic (ROC) curve analysis, in order to select the rank aggregation that shows better concordance with literature-derived toxicity genes (see Methods, subsection ROC analysis for details). The optimal rank aggregation was used to generate four gene candidate lists, one for each model system, which were finally subjected to functional enrichment analysis to evaluate the ability of the four models (HH, RH, RL, RK) to recapitulate human organ-specific drug toxicity.

Differential gene expression analysis

The overall analysis across the four model systems shows that the literature-derived toxicity-related human genes are more often dys-regulated in the in vitro scenarios as compared to the in vivo ones (Fig 1). On average, we have far more genes dys-regulated in both in vitro systems (HH and RH) as compared to both in vivo systems (RL and RK), which might be a consequence of isolation stress or because of the cell culture conditions as noted before by Godoy et al. [39].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Summary of differential expression analysis.

Distribution of the number of differentially expressed genes is depicted with violin-plots by model systems, with each dot representing the number of genes modulated by the high dose treatment with a single compound. The median and interquartile range of each distribution is depicted by black squares and black bars, respectively. On the same plot, we illustrate the distribution of the full set (all), the filtered set with unique human-rat orthologs (ortholog) and the subset found in the corresponding organ-specific literature-derived genes (literature).

https://doi.org/10.1371/journal.pone.0210467.g001

Furthermore, we observe a large discrepancy in the number of differentially expressed genes across compounds, which shows good correlation with the increase of dose level (S1 Fig), and in many cases, also with the known toxicity label of the drug. For example, indomethacin is a known toxicant with high liver injury risk and shows high number of differentially expressed genes in the hepatocyte and liver systems (S2 Fig). The same holds for the hepatotoxic analgesic acetaminophen (commonly known as paracetamol), although we observe an inconsistent pattern as well, e.g., known hepatotoxic drugs such as enalapril and valproic acid dys-regulated ten (or more) times less genes in the RL than in the RH system.

Rank aggregation across compounds and scenarios

The ROC evaluation of the aggregation methods (Fig 2A, S3 Fig) across compounds suggests that the rankings from models with largest numbers of differentially expressed genes recapitulate better the genes from the literature. Overall, in the in vitro scenario we get better performance (larger AUC) as compared to the in vivo one. While the highly imbalanced size of differentially expressed gene sets seems to affect the MEAN aggregation method the most, the NC rank, which simply uses the number of compounds modulating a given gene, works in par with the more complex aggregation methods. Due to the influence of the partial ranking, the MEAN method and to a lesser extent the RRA method delivered often the same ranks for the larger sets of genes; at the same time, the NC rank was quite high for these. The latter indicated that we might benefit if we do second-level aggregation across the rank aggregation methods, and one simple way is to go for the median of the aggregated ranks (named AR).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Performance comparison of gene ranking methods.

A) ROC curves for compound-wise rankings (gray curves) and compound-aggregated rankings (coloured curves) as calculated by the rank-aggregation methods, with literature-derived genes used as positive examples. Curves are not smoothed, and the straight line fragments of the curves result from the initial partial ranking (same lowest rank for all not-modulated genes, see Methods for details). B) Summary of AUC (area under the ROC curve) values for all rankings across all model systems, including model system-wise median AUC represented as dashed line. AUC values above 50% means better performance than randomly assigned rank.

https://doi.org/10.1371/journal.pone.0210467.g002

Evaluation of aggregation methods indicated that STUART was the optimal with AR very close (Fig 2B). Since AR is an aggregation method, we preferred it for the final gene prioritization. We selected the top 400 AR-ranked genes in each scenario for down-stream analysis.

The four rankings reflect how the four systems (HH, RH, RL, RK) respond to the same set of toxic compounds and thus allow us to evaluate their strengths and their limitations in relation to each other in terms of recapitulating different aspects of human organ toxicity as reflected in affected gene expression. The connection to genes allows the generation of hypotheses about the concrete biological processes involved.

Examination of the highest 10 ranked genes (Table 1, columns HH, RH, RL, RK) shows no overlap between models. A more meaningful result is obtained when fully comparing the four rankings (Fig 3): the best correlation found is between RH and RL (Spearman’s rho = 0.61). RH to HH is second best (Spearman’s rho = 0.49), which suggest that they certainly differ in terms of response to toxic drugs, but not so much than when comparing other models (e.g., RK versus RL, Spearman’s rho = 0.28).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Comparison of AR-ranked gene sets across model systems.

Rank correlation between the model system-specific gene rankings as measured by the Spearman’s rho correlation coefficient.

https://doi.org/10.1371/journal.pone.0210467.g003

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Top 10 ranking genes for each model (HH, RH, RL and RK) and for combined rankings: Sum (all), and subtractions (HH>RH, RH>HH, RL>RH and RH>RL).

Font indicates toxicity association by literature, bold for liver, underlined for kidney.

https://doi.org/10.1371/journal.pone.0210467.t001

To find out about the processes affected by toxicity in each of the four models, we did functional enrichment analysis of various terms for which the genes were annotated (including Gene Ontology terms, pathways, participation in protein complexes, and tissue expression; see Methods for details).

Interestingly, now the similarity between the top-100 genes of HH and RH is manifested in terms of functional GO BP (Gene Ontology Biological Process) terms, where the most significant are “carboxylic acid metabolic process” and “monocarboxylic acid metabolic process”, respectively (see S5 and S6 Files). In both systems “response to drug” appears, although at a less significant level.

A similar analysis for the top-100 genes of RL indicates terms related to cell cycle and cell division (S7 File) indicating that the processes most affected are less specific when studying an organ. In the case of RK we did not even obtain terms at high significance level (P-value < 1e-03; S8 File).

Discussion

Here, we presented a protocol to extract ranked gene lists from toxicogenomics datasets specific to model systems (e.g. human hepatocytes), which includes evaluation by agreement with the literature in human toxicity (e.g. dealing with human hepatotoxic genes). This approach allows evaluating different ranking methods by the agreement of the resulting gene rankings to the literature, and can be used to evaluate and compare the aspects of human toxicity recapitulated by each system model.

Due to the multi-factorial experimental designs of the existing toxicological databases, we needed to integrate data across multiple model systems, drugs, dose levels and time points. A relatively simple way is to summarize the genome-wide drug-induced transcriptional changes in a prioritized (ranked) gene list, followed by robust rank-aggregation approach [15] to combine them into a single organ-specific gene list. While this approach is commonly used for prioritization of disease-related genes [14], here we used it as baseline method for prioritization of toxicity-related genes. In particular, we exploit the robustness (with respect to noise, outliers, errors) and computational efficiency of recently proposed probabilistic models for rank-aggregation, i.e., RRA [15] and BIRRA [16].

While the task of aggregating results from toxicogenomics datasets is clearly not trivial [3], current research efforts show that carefully designed protocols for meta-analysis can potentially address issues of heterogeneity of samples. For example, novel (more generic) gene biomarkers for drug-induced toxicity have been identified by integration of differential expression analysis and data mining of in vitro and/or in vivo toxicogenomics data [47, 48]. Several studies exploited biclustering methods [49] for coherent coexpression gene module discovery combined with downstream functional enrichment analysis to characterize conserved pattern of drug response across human cell lines [50] and prioritize gene modules for specific chemically-induced liver [51] and kidney injury endpoints [52]. Last, but not least, El-Hachem et al. [53] showed that mapping the drug-induced transcriptional changes in hepatocytes into the space of enriched biochemical pathways, combined with subsequent bi-clustering, reveals conserved modules of chemically modulated pathways in human hepatocytes, rat hepatocytes and rat liver model systems.

Here, we focused on the evaluation of organ-specific toxicity-associated genes, and proposed a protocol to prioritize genes that we evaluated with the literature data, with a focus on kidney and liver. We applied this approach to a toxicogenomics dataset corresponding to 33 compounds (mainly toxic drugs) profiled in rat kidney (RK), rat liver (RL), primary rat and human hepatocytes (RH and HH, respectively), obtained from Open TG-GATEs. This resource was especially fitting our purpose, due to consistent experimental design across the four target systems; TG-GATEs has been recently exploited to find hepatotoxic compounds [54] and to compare rat liver testing systems [55]. We evaluated multiple ways to aggregate these data for the description of drug-induced toxicity in human kidney and liver.

In our analysis, we chose maSigPro for modeling gene expression as it perfectly fits the dynamic nature of the data (multiple time series). On the other hand, maSigPro selects significantly different gene expression profiles based on two-step approach, in which genes that do not pass the first step (for a given significance) will not be modeled in the next step (no P-value and R2 calculated), meaning a worst possible rank in our case. This means that the compound-wise gene lists are technically partially ranked lists. While this can be a problem for standard rank-aggregation methods (as evident for the mean rank in the ROC analysis), the alternative approaches RRA, Stuart and BIRRA were not specifically affected, making the proposed method suitable for meta-analysis of partially ranked gene lists across compounds obtained from literature.

There is a substantial level of gene filtering in our protocol, which was necessary to fairly compare the rank-aggregation methods across species and with the literature-derived human gene lists. The current method uses only one-to-one orthologs as defined in NCBI Homologene databases. Furthermore, the current approach for generating toxicity-related genes from annotated biomedical literature takes only into account human protein-coding genes. This means that the results should improve by inclusion of more orthology relations (available in other databases) as well as by extension of the literature-derived lists with toxicity related non-coding human genes.

The differential expression analysis as well as the rank-aggregation analysis across drugs show that in the in vitro model systems we get better coverage of the human gene candidates from the literature as compared to the in vivo one. The highly imbalanced size of the sets of differentially expressed genes seems to affect the mean aggregation method the most. Despite this, aggregation of partial rankings is still capable to find common drug-response genes across all model systems. An additional factor that could explain the better agreement to the in vitro models is the fact that toxicity is often tested in cell lines and not in vivo. For example, around a 30% of references related to the toxicity genes ATF3, GDF15, HMOX1 and EGR1 mentions also a cell line.

Another factor that could influence the differences between the ranks derived from the four model systems for their agreement with the literature-derived lists could lie in the differences in numbers of genes associated in the literature (which itself could reflect real physiological differences). Accordingly, the lowest agreement found in RK could be related to the fact that the number of literature-derived genes associated with kidney toxicity was much smaller than those associated to liver toxicity. Both the HH and RH system recapitulated better the literature than RL, maybe indicating that the in vitro systems are more appropriate in general to model liver toxicity. Unexpectedly, RH was a bit better than HH: while the similar performance indicates that RH is a good model for studying liver toxicity, this result could be influenced by the fact that studies dealing with human genes (among the selected PubMed abstracts used for extraction of the literature-derived gene lists) might actually use rodent models.

Ranked lists for the four models, in vitro HH and RH, and in vivo RL and RK, indicated little agreement in the top ranked genes (Table 1). Full rank comparisons (Fig 3) suggested that the two rat hepatic models, RH and RL, are the most similar, while the two in vitro hepatocyte models (HH and RH) follow next. Functional analysis, however, suggested that HH and RH were the models more similar in terms of the functions of genes affected, mostly related to carboxylic acid metabolism.

Next, we combined the ranks looking for genes ranking high in the four models. Adding pragmatically all ranks brings HMOX1, a very well-known toxicity associated gene to the top position, and highlights a significant number of genes related to responses to toxicity and chemical stimuli, suggesting that the addition of rankings is meaningful. Collectively, cell cycle and cell division seem to be the processes most generally affected in the models. Rank comparison between HH and RH highlights genes with astonishing different reactions to toxicity. This suggest that researchers should not use RH as toxicity model if they are interested in certain genes related to human liver toxicity that react in HH but are non-responders in RH, such as GSTP1. Similarly, our comparison of RH and RL indicates that while researchers interested on toxic effects of genes related to lipid binding should be using hepatocytes, important toxic effects on transporters that are observed in RL will be missed if using RH.

As every integrative approach that relies on multiple subsequent steps of data preprocessing, modeling, and high-level summarization, our integration protocol can be influenced by the limitations of each subsequent step of the analysis: from comprehensiveness of the input (toxicogenomics) data, issues with data normalization and gene expression modeling, to specificity and sensitiveness of the method (maSigPro) to detected gene expression changes, the criteria used for ranking (P-value and R2 of the gene-wise models) as well as the robustness of the subsequent aggregation of gene lists (median of ranks, rank sum). However, the proposed protocol is fully flexible in the sense that each of the subsequent steps can be easily replaced with suitable alternatives.

Our approach to find literature derived genes also has some limitations. For example, we missed PON1, which is known to be related to liver toxicity (see e.g., [56]). The reason is that we used a conservative approach with stringent cut-offs in the selection of positives (see Methods for details).

Here we have presented a protocol for the aggregation and literature evaluation of toxicogenomics data. Application to liver and kidney model systems of in vitro and in vivo toxicity allowed us to define genes and pathways with general and particular toxicity related roles in these models. We believe that this protocol should be easily applicable to other human organs and to other models in toxicogenomics datasets.

Supporting information