Simultaneous Non-Negative Matrix Factorization for Multiple Large Scale Gene Expression Datasets in Toxicology

Non-negative matrix factorization is a useful tool for reducing the dimension of large datasets. This work considers simultaneous non-negative matrix factorization of multiple sources of data. In particular, we perform the first study that involves more than two datasets. We discuss the algorithmic issues required to convert the approach into a practical computational tool and apply the technique to new gene expression data quantifying the molecular changes in four tissue types due to different dosages of an experimental panPPAR agonist in mouse. This study is of interest in toxicology because, whilst PPARs form potential therapeutic targets for diabetes, it is known that they can induce serious side-effects. Our results show that the practical simultaneous non-negative matrix factorization developed here can add value to the data analysis. In particular, we find that factorizing the data as a single object allows us to distinguish between the four tissue types, but does not correctly reproduce the known dosage level groups. Applying our new approach, which treats the four tissue types as providing distinct, but related, datasets, we find that the dosage level groups are respected. The new algorithm then provides separate gene list orderings that can be studied for each tissue type, and compared with the ordering arising from the single factorization. We find that many of our conclusions can be corroborated with known biological behaviour, and others offer new insights into the toxicological effects. Overall, the algorithm shows promise for early detection of toxicity in the drug discovery process.


Introduction
The aim of this work is to highlight the usefulness of a recently proposed extension to the technique of non-negative matrix factorization (NMF) by demonstrating its promise for early detection of toxicity in the drug discovery process. In particular, we (a) show that any number of related datasets can be treated simultaneously with this approach, (b) deal with practical issues that arise when the algorithm is applied to real datasets, (c) demonstrate its use with a new large scale microrray dataset, and (d) interpret the results from a biological perspective.

Computational Background
NMF seeks to represent a large complex dataset in terms of smaller factors. The name covers many algorithms. Each approximates a non-negative matrix as the product of two or more smaller non-negative matrices, by attempting to minimise some objective function. Lee and Seung [1] showed that when applying multiplicative non-negative factorization to images of faces, each row/column pair of the factors expresses a recognisable facial feature. These techniques have since been used in many settings to learn parts of the data as well as to factorize and cluster datasets. For example, when applied to text data in [1] the algorithm can differentiate multiple meanings of the same word by context. On microarray data, NMF has been used to find patterns in genes or samples, typically bi-clustering both groups in a similar manner to two-way hierachical clustering [2][3][4][5][6][7]. The review article [8] shows how NMF has also been successful in other areas of computational biology, including molecular pattern discovery, class comparison and biomedical informatics. The new challenge that we address in this work is to apply the NMF methodology to multiple, related, large scale, data sets simultaneously. We use the The area under consensus cumulative density, [3,14]. (c) The cophenetic correlation coefficient, [3]. doi:10.1371/journal.pone.0048238.g001 Table 1. Blood clinical chemistry analysis for each mouse.  Vehicle  42  188  12  4  348  484   I  B  Vehicle  41  92  9  6  364  258   I  C  Vehicle  29  75  9  6  278  166   II  D  6  95  441  11  8  1218  4930   II  E  6  692  981  8  7  2126  1130   II  F  6  52  83  9  8  294  152   III  G  20  312  1300  6  8  3172  2544   III  H  20  462  937  8  6  1760  1182   III  I  20  698  1090  6  7 2616 1592 The mice were randomly divided into three groups and treated with either Vehicle or two concentrations of PPM201 (6 or 20 mg/kg body weight). The response to the ''therapeutic dose'', 6 mg/kg, was found to vary widely for ALT (alanine aminotransferase), AST (aspartate aminotransferase), LDH (lactate dehydrogenase) and CK (creatine kinase). AST is raised in PPM201 treated animals, with mouse E (6 mg/kg) seeming to be especially raised; AST is known to be variable between animals, but mouse E also shows a higher level of ALT, indicating that there may be a shared mechanism for the two enzymes. Creatinine is decreased in liver and possibly kidney disease; the contrasts observed here are inconclusive. BUN (Blood, Urea and Nitrogen) is raised in kidney disease; results are again inconclusive. Following cardiac infarction LDH is increased after 12 hours, possibly also caused by liver toxicity; mouse E is markedly lower than the other PPM201 treated animals and it may be that its heart muscle profile might be more similar to the untreated mice. CK is, like LDH, increased in myocardial infarction and this supports the LDH findings for mouse E. doi:10.1371/journal.pone.0048238.t001 Figure 2. Factorising as a single dataset; reordering using the NMF for k~4. The columns show the samples and the rows the gene expression for each of the 45037 genes. Genes and samples are organised by cluster number. Elements within each cluster are ordered, with the largest value at the bottom/right. Each tissue is characterised by a group of highly expressed genes; from the top left to bottom right these are heart, skeletal muscle, liver and kidney. For comparison purposes, the characteristic 100 ''best'' genes in the four columns are names heart 1 , skeletal muscle 1 , liver 1 and kidney 1 . doi:10.1371/journal.pone.0048238.g002 work of Badea [9,10], who considered an extension of NMF that deals with two data matrices. Simultaneous NMF is used in [9] to study pancreatic cancer microarray data alongside extra information concerning transcription regulatory factors. In [10] microarray datasets for pancreatic ductal adenocarcinoma and sporadic colon adenocarcinoma are sumiltaneously factorized in order to discover expression patterns common to both data sets. This simultaneous NMF approach readily extends to the case of an arbitrary number of data matrices and here, for what we believe to be the first time, we implement and evaluate the method on more than two. We also consider various practical issues that must be tackled in order to produce a useful computational tool. To minimize the number of algorithmic parameters, make the results straightforward to interpret, and exploit the natural sparsity in the algorithm [9, section 3], we focus on hard clustering. The interesting issue of allowing clusters to overlap in this context is therefore left as future work.

Biological Background
We analyse gene expression data describing the molecular changes in four tissue types due to different dosages of an experimental pan-peroxisome proliferator-activated receptor (pan-PPAR) agonist PPM-201, provided by Plexxikon. PPARs have attracted great interest as potential therapeutic targets for diabetes [11], but major concerns have arisen due to clinically observed side-effects [12]. Hence, there are compelling reasons for toxicological studies at the gene expression level.
The material is organised as follows. In Section we describe the simultaneous NMF algorithm and outline our approach for using the output to order and cluster a dataset. Section describes the mouse microarray data, and the NMF results that arise when we treat it as a single dataset are given in Section . This is followed in Section by the analysis of the data split into four datasets corresponding to the known tissue types; liver, kidney, heart and skeletal muscle. In Section we compare the gene clusters from Sections and , and Section discusses the results. Conclusions are given in Section .

Algorithms
Given d non-negative data matrices A(i) of size m(i)|n for i~1, . . . ,d, our aim is to simultaneously factorize all matrices so that with the additional constraints that W (i) is a non-negative matrix of size m(i)|k for i~1, . . . ,d, and H is a non-negative matrix of size k|n. Generalising naturally from the d~2 case in [9], we seek to minimise the objective function Here E : E denotes the Frobenius norm. As in [9] the b coefficients are designed to give equal weight to the different error terms. Based on the multiplicative update rules developed in [13], an iterative algorithm that attempts to solve the optimisation problem can be derived using a gradient descent method dz1 times. This gives us the following sequence of approximations for j~1,2, . . ., given initial choices W (i) ½0 and H ½0 , for some small positive matrices g ½j W (i) , and g ½j H , with Ã representing element-wise multiplication. The iteration may be motivated through the intuition that when g ½j W (i) and g ½j H are sufficiently small and positive each of these equations should reduce the objective function. This allows us to set again with the division being performed element-wise. Hence the overall iteration has the form The values in g W (i) and g H are non-negative due to the constraints on the matrices, however they are not necessarily small. The iteration decreases the objective function (1), so this leads to a locally optimum solution, but we cannot guarantee convergence to a global optimum. In particular, different initial conditions can lead to different factorizations of different quality.
Having iterated up to some stopping criterion and produced the factorizations, we use them to bi-cluster the data. Each sample is assigned to the cluster for which it has the largest value in the gene cluster and vice versa. In reordering the data for easy visualisation we organise the rows and columns by cluster number (assigned arbitrarily) and sort the elements within each cluster from the appropriate sample/gene set, with the largest value at the bottom/ right of that cluster. Given that the second factor is common to all the factorizations, it produces a matching ordering of the columns of the data.
Because the result depends on the choice of initial condition, and because the choice of k is not automatic, further information is needed in order to specify a practical algorithm. To deal with the lack of uniqueness, we try several initial conditions and pick a realisation that minimises the objective function (1). We then continue until further runs do not significantly alter the results. The objective function value is also one of the criteria we use in order to decide which rank/clustering is the most ''appropriate'' for the data. By regarding the objective function as a function of k, we identify values of k where the decay in the objective function begins to diminish. In addition we also form a consensus matrix as in [3,14] for the clustering of the objects. This is the average of the connectivity matrices C where for each initialisation C i,j~1 if objects i and j are clustered together and 0 otherwise. So the consensus matrix contains values between 0 and 1 with the (i,j)th element being the likelihood that objects i and j cluster together. The cumulative density of these values is constructed, by summing the appropriate probabilities, and the area under this curve is the second measure we look at when considering choices for k. The third measure is the Pearson correlation of the cophenetic distances, as explained in [3].

Mouse data
We apply these techniques to mouse gene expression data quantifying changes in four different tissue types following administration of different dosages (vehicle, therapeutic and toxic) of an experimental pan-PPAR agonist. The study design and clinical chemistry results are summarised in Table 1. ALT and AST are known markers in rodents for liver toxicity [15] and from this criterion mouse E may be showing a toxic response to The area under consensus cumulative density function for k~2, . . . ,10, [3,14]. (c) The cophenetic correlation coefficient, [3]. doi:10.1371/journal.pone.0048238.g004 PPM201, despite it being administered at a supposedly therapeutic dose level. This conditions our expectation of the gene-expression pattern for mouse E and suggests that it may be similar to the toxic level group III for liver.
Nine wild type mice (strain: C57BL/6J) were randomly divided into three groups; -Group-I, II and III. PPM-201 in the vehicle base was administered daily for 14 days at 6 mg/kg body weight dose rate to each mouse in Group-II and at 20 mg/kg body weight dose rate to each mouse in Group-III while the mice in Group-I received only the vehicle base. On 15th day, the mice were sacrificed to harvest blood, heart, skeletal muscle, liver and kidney tissues for clinical chemistry, microarray and histopathology Figure 5. Factorisation of the four separate tissue types using simultaneous NMF with k~3. Top left, kidney; top right, liver; lower left, heart; lower right, skeletal muscle. The four tissue types are treated as separate sources of information across a common set of mice. Genes are therefore ordered differently in each of the four tissues, but the mice ordering is global. The resulting mouse ordering and mouse clusters are detailed in Table 3. doi:10.1371/journal.pone.0048238.g005 analysis. In the clinical chemistry analysis, alanine aminotransferase (ALT, U/L), aspartate aminotransferase (AST, U/L), creatinine kinase (CK, U/L), blood urea nitrogen (BUN, mmol/L), creatinine (mmol/L) and lactate dehydrogenase (LDH, U/L) were measured from the blood of each mouse. Two sections of liver, two sections of kidney, one or two sections of skeletal muscle, and one section of heart were prepared from each mouse, stained with hematoxylin and eosin (H&E), and examined by a veterinary pathologist. Total RNA was isolated from murine tissues using Qiazol-based homogenization and subsequent column-based purification (Qiagen) with on-column DNAse-treatment. DNAsefree RNA was assessed for quality using Agilent Bioanalyser electrophoresis and acceptance criteria of RNA Integrity Number (RIN) greater than seven. 50 ng of total RNA was subsequently utilized as input to cDNA-based amplification and biotin-labelling using single-primer isothermal amplification according to the manufacturer's instructions (Ovation System, NuGEN Technologies). Unlabelled and biotin-labelled cDNA was qualitatively assessed by Agilent Bioanalyser electrophoresis to ensure identical size distributions of all samples pre-and post-fragmentation. Fragmented, biotin-labelled cDNA were hybridized to MOE430 2.0 GeneChip arrays (Affymetrix) with subsequent scanning and feature extraction according to the manufacturer's instructions.
The dataset has been approved by the GEO curators and assigned the accession number GSE31561.   Table 3. Ordering of the tissue samples after a four-way factorization of rank 3.

Single dataset
First, the samples are treated as a single dataset, with thirty six samples and 45037 genes, hence the data matrix A is 45037|36. This corresponds to the case where d~1 in Section . The factorizations were performed twenty times for each k~2,3, . . . ,16, with a consensus matrix formed from the clustering of the samples. All gene clusters associated with this analysis are labelled with a subscript 1, e.g., heart 1 . Figure 1(a) shows the minimum size of the objective function that we observed for each value of k. We see that this value decreases monotonically, with a slower rate starting at around k~4. Figure 1(b) shows the area under the cumulative density curves for the same values of k. This subfigure clearly points to k~4, as does subfigure (c) showing the cophenetic correlation.
Based on Figure 1, we conclude that when the data is factorized as a single entity, k~4 clusters is the most appropriate choice. Reordering the dataset using the ordering for k~4 in the manner described in Section gives the images shown in Figure 2. This figure shows the samples in the columns with cluster one at the top. To aid visualisation, the sample clusters are split by white lines, as are the gene clusters. This reordered data matrix shows a distinctive ''ramp'' effect in the blocks on the diagonal, placing genes that are most influential in identifying each tissue type to the bottom of the block. This figure also shows some of the differences in expression behaviour between the tissue types, particularly for the most influential genes.
Because we know the origin of the samples, we can confirm that the algorithm has put the heart samples in cluster one, the skeletal muscle samples in cluster two, the liver samples in cluster three, and the kidney samples in cluster four. The exact ordering of the samples is shown in Table 2. This table also shows the mouse identification information for each sample, and we see that the mice are not ordered in the same way within each cluster. It is the liver and skeletal muscle samples that most closely respect the dosage levels within the clusters. Both these clusters only have one sample mis-ordered.
Given that the factorization has been performed for k~2, . . . , 16 we know what the clustering would be from all these rank factorizations. This information is displayed in Figure 3. Here the rows representing the samples are ordered in tissue then dosage subgroups. For each rank k, samples with the same colour are assigned to the same cluster. As we have seen before, for k~4 the samples are split into tissue types. The figure shows that this split persists at k~5 with an empty cluster forming. In fact, for this range of k there are at most twelve clusters of samples. We also see from this figure that for no value of k are the twelve tissue/dosage subgroups found.

Multiple datasets
The test in Section indicates that the basic NMF factorization approach can deliver biologically meaningful results-separating the twelve samples by tissue type. But the failure to order correctly within tissue type according to dosage motivates the use of the multiple dataset generalization introduced in Section , where the four tissue types are treated as separate sources of information across a common set of mice. Intuitively, we would expect to add value to the data analysis by building known biology into the algorithm in this way. In this section, we therefore factorize the four new datasets simultaneously. This is similar to the test in Section in the sense that it produces a single ordering for the mice, but it has the potential to add extra information by providing four different, tissue-level, gene orderings. We thus have d~4 matrices H1, SM1, L1 and K1 are the gene clusters most characteristic for the heart, skeletal muscle, liver and kidney, respectively, in the single (combined) data set, as in Figure 2.
Clust.1, 2, or 3 denotes the 100 genes most securely placed within the clusters of the diferently ordered genes in the 4-way factorization shown in Figure 5. The order of the clusters is 1-3, from the top of the figire, for each tissue. We refer to these clusters as ''heart 4 cluster1,'' etc. The overlap of the heart 1 from the one-way factorization to heart 4 is referred to as heart 1 heart 4 cluster 1. doi:10.1371/journal.pone.0048238.t004 of size 45037|9. We again performed 20 factorizations, this time for k~2, . . . ,10 and these have been used to generate a consensus for clustering the mice. The objective function and the consensus measurements are shown in Figure 4. The objective function in subfigure (a) does not show much decrease in convergence rate until we get to nine clusters. This is the point where each mouse is put into a cluster on its own. The area under the cumulative density curve in Figure 4(b) suggests using either rank k~3, or k~5 factorizations for the clustering. The correlation coefficients shown in subfigure (c) give the same two values as peaks, as well as k~8, though the k~3 peak is the highest.
Given these measurements we consider the four-way simultaneous factorization for k~3 in Figure 5. The reordered datasets are shown separately with the kidney dataset in the top left, the liver dataset in the top right, the heart dataset in the bottom left and the skeletal muscle in the bottom right. The mouse ordering and mouse clusters that arise are shown in Table 3. The four subfigures in Figure 5 also illustrate that the gene clusters are different for each dataset. The three clusters for each tissue in this 4-way factorization are subsequently refered to in the form ''tissue 4 , cluster 1,2 or 3. '' Table 3 shows that the simultaneous NMF approach has recovered the known mouse treatments except for one misplacement. Figure 6 shows the clustering for the fourway simultaneous factorizations for k~1, . . . ,9. This indicates that this mouse does not cluster with all those of the same dosage for any rank of factorization greater than two, instead it associates with the higher more toxic dosage. This is borne out by the known blood chemistry, as summarised in Table 1; the mouse that is mis- Figure 7. Enrichment of canonical pathways in the four tissue specific gene clusters. The top one hundred most influential probe-sets in the four tissue specific gene clusters obtained in the first factorization were subjected to signalling and metabolic pathways analysis in the IPA software. This graph shows the comparison of canonical pathways enriched in the four tissue specific gene clusters, heart 1 , muscle 1 , kidney 1 and liver 1 . The coloured bars show the significance of the enrichment for a particular pathway in the cluster computed by Fisher's exact test. doi:10.1371/journal.pone.0048238.g007 classified exhibits a toxic response and is therefore classified with the mice that received the higher dose.

Comparing Gene clusters
Our aim now is to test the results from the novel multi-way NMF algorithm used in Section in order to see whether they (a) show consistency and (b) add value to the results in Section from standard NMF. We know that the four simultaneously factorized datasets correspond to the four clusters of samples that were discovered in an unsupervised manner from the single factorization of the full dataset. It could therefore be conjectured that the most influential genes in the first factorization will appear as influential genes in the four-way simultaneous factorization for that dataset, but less so for the other datasets.
Our comparisons involve four reference sets. For illustration, we chose an arbitrary threshold of one hundred; that is, we consider the top one hundred most influential genes from the four clusters in the first factorization shown in Figure 2. For easy reference these sets are referred to using the known tissue type. This means that the genes from cluster one are the heart 1 genes, those from cluster two are the muscle 1 genes, those from cluster three are the liver 1 genes and those from cluster four are the kidney 1 genes. The 4-way factorization shown in Figure 5 identifies differently ordered gene clusters for each tissue, which we will refer to as ''kidney 4 , cluster 1,2 or 3, etc. '' Table 4 shows the total number of co-incident genes between the top 100 lists arising from the oneway and four-way factorisations. The table also shows the probability of the two lists having that number of genes in common if the second list were randomly selected; hence these values come from the hypergeometric distribution. We see that the important genes for each tissue type appear significantly highly in the clusters from that tissue's data type. In addition, all the tissue type genes also appear significantly within the reordering of the heart dataset. This link is reciprocated, with the heart genes appearing significantly frequently within the skeletal muscle dataset. Surprisingly, the greatest overlap arose between liver 1 and heart 4 cluster 2. One of these genes, Apoliprotein A1, is being considered as a marker for cardiac toxicity [16].
We would like to demonstrate the utility of the factorization method by using the gene clusters obtained in our analysis to understand tissue specific effects of the experimental drug, PPM-201. Of course, we are not claiming that this is an exhaustive analysis of the effects of PPM-201. We analysed the gene clusters for pathways enrichment and gene ontology enrichment using DAVID [17] and Ingenuity Pathways Analysis (IPA) [18] tools. Table 5 shows the comparison of KEGG pathways enriched in the four tissue specific top one hundred most influential probe-sets obtained in the first factorization. Pathways enriched in these clusters differ according to the tissue types and can be considered as the pathways that are most perturbed by PPM-201. For example, arrhythmogenic right ventricular cardiomyopathy, hypertrophic cardiomyopathy and dilated cardiomyopathy are enriched in heart, whereas starch and sucrose metabolism, drug metabolism and PPAR signalling pathway are enriched in liver. Similarly, Figure 7 shows the enrichment of canonical pathways in the four tissue specific clusters analysed using IPA. It also shows the tissue specific enrichment of pathways-calcium signalling, integrin linked kinase (ILK) signalling and cardiac hypertrophy signalling are enriched in heart 1 and muscle 1 clusters, whereas fatty acid metabolism and farnesoid X receptor (FXR)/retinoid X receptor (RXR) activation are enriched in the liver 1 cluster.  Table 6. Enrichment of KEGG pathways in the common genes between the clusters found by the two ways of factorization.    Analysis of the same sets of genes for enrichment of toxicity functions in the IPA shows, in Figure 8, cardiac hypertrophy in heart 1 genes, increased level of creatinine and hydronephrosis in kidney 1 genes, and increased levels of lactate dehydrogenase (LDH) and steatosis in liver 1 genes.
The common genes between the top one hundred most influential probe-sets in the four tissue specific clusters and the top one hundred most influential probe-sets in the clusters formed by 4-way simultaneous factorization of the split dataset were also analysed for enrichment of pathways, gene ontology and toxicity functions using DAVID and IPA. Tables 6,7,8,9,10,11,12,13 ,   Table 10. Heart 1 heart 4 cluster 1. Common probesets between the top one hundred most influential probesets in the heart 1 cluster and 20 mg/kg dosage cluster (cluster 1) of the heart 4 dataset.

Sr.
Probeset ID Gene Symbol Entrez Gene Name 14, 15, 16, 17 summarise the results of this analysis, which are discussed further in the next section.

Discussion
The factorization and reordering of the dataset as a whole set ( Figure 2 and Table 2) successfully clustered samples from the same tissue and further investigation showed that it simultaneously identified genes with a known relevance to those tissues. It was therefore reasonable to study the genes that were responsible for this differentiation. In the one-way clustering, the top 100 probesets from each of the four tissue specific clusters show remarkable coherence for tissue specific pathways. The calcium signalling pathway is highly enriched in both heart 1 and muscle 1 clusters; these genes are linked to muscle contraction function. Muscle contraction is the prime function of cardiac and skeletal muscles. A deeper look at the probe-sets (Tables 7 and 8) from the heart and skeletal muscle clusters shows a successful identification of differences in the tissue types for this pathway; see Figure 9. MYH1, MYH2, MYH4 and MYL1 of the myosin family, which are specific to skeletal muscle, are found in the muscle 1 cluster while cardiac muscle specific myosin family members MYH6, MYL2 and MYL3 are found in the heart 1 cluster [19]. This pattern is also true for troponin, calsequestrin, ryanodine and actin     of muscle 1 cluster and 20 mg/kg dosage cluster of muscle 4 , and between the top 100 probe-sets of muscle 1 and vehicle dose cluster of muscle 4 , 15 and 14 probe-sets were in common and are named as muscle 1 muscle 4 cluster 1 and muscle 1 muscle 4 cluster 3, respectively. The calcium signalling-skeletal muscle contraction pathway is enriched in muscle 1 muscle 4 cluster 1 with the presence of ACTA1, ATP2A1, MYH1, MYH4, RYR1, TNNC2, TNNI2, TNNT3 and TRDN genes, whereas muscle 1 muscle 4 cluster 3 does not show any significant enrichment for signalling or metabolic pathways. Interestingly, 49 probe-sets in the liver 1 heart 4 cluster 2 are common between the top 100 probe-sets of liver 1 cluster cluster and 6 mg/kg dosage cluster of liver 4 and highly enriched for acute phase response signalling, prothrombin activation and FXR/RXR activation pathways with the presence of ALB, ABCB11, AMBP, APOA1, APOA2, APOC3, APOH, F2, FGA, FGB, FGG, HAMP, HP, HPX, ORM1, PON1, RBP4, SERPINA1, SER-PINC1, SLC27A5 and SLCO1B3 genes ( Figure 11). This suggests alterations in lipid metabolism in liver along with tissue injury in heart induced by PPM-201 at 6 mg/kg dosage [30][31][32][33], which becomes more plausible when we look at the genes in liver 1 liver 4 cluster 1 that are common between the top liver 4 genes and 20 mg/kg dosage cluster of liver 4 dataset. Enrichment of toxicity functions in liver 1 heart 4 cluster 2 using IPA shows increased level of LDH as one of the toxicity functions ( Figure 12) which has been validated with the increased level of LDH in the clinical chemistry results.

Conclusions
We have demonstrated that multi-way simultaneous nonnegative matrix factorization can be usefully applied to the case of multiple datasets-here, for what we believe to be the first time, more than two large scale matrices were treated. The results were shown to be consistent with, and to add value to, standard nonnegative matrix factorization of the whole dataset.
In summarizing our biological findings, we first note that the roles of the three different isoforms of PPARs -PPAR-a, PPAR-b (also known as PPAR-d) and PPAR-c in metabolism and their difference in expression in different tissues and different species are well known [34][35][36]. In mouse, PPAR-a is highly expressed in liver and to a lesser degree in kidney, heart and skeletal muscle; PPARb is expressed in many tissues but peaks in kidney, heart and intestine whereas PPAR-c is mostly expressed in adipose tissue Figure 10. Liver genes enriched in FXR/RXR activation pathway IPA analysis of the top 100 probe-sets from the liver 1 cluster (Figure 7) showed the enrichment of FXR/RXR activation pathway. The genes present in this pathway are highlighted in orange. The liver genes present in the pathway are given in Table 9. Pathway diagram was drawn using Path Designer function of IPA [18]. doi:10.1371/journal.pone.0048238.g010 Figure 11. Enrichment of canonical pathways in the liver heart gene cluster no. 2. This gene cluster has 49 common probe-sets between the top one hundred most influential probe-sets in the liver gene cluster and top one hundred probe-sets in cluster number 2 (6 mg/kg dose rate) of the heart dataset reordered by 4-way simultaneous factorization. Canonical pathways enrichment for these 49 probe-sets analysed using the IPA software is shown in this figure. The length of the bars shows the Fisher's exact test p-value for enrichment for a particular pathway in the cluster. doi:10.1371/journal.pone.0048238.g011 Figure 12. Enrichment of toxicity functions in liver 1 heart 4 cluster 2. This gene cluster has 49 common probe-sets between the top one hundred most influential probe-sets in liver 1 heart 4 cluster 2 (6 mg/kg dose rate). Toxicity functions enrichment for these 49 probe-sets analysed using the IPA software is shown in this figure. The length of the bars shows the Fisher's exact test p-value for enrichment for a particular pathway in the cluster. doi:10.1371/journal.pone.0048238.g012 [34,37]. Pan-PPAR agonists activate two or all of the pan-PPAR isoforms and differ in their pharmacological actions. Factorisation of the dataset after splitting it on the tissue basis appears to be beneficial in identifying tissue specific and dosage effects of the experimental pan-PPAR agonist PPM-201 in this study. This approach could be useful in understanding molecular mechanisms and identifying potential tissue specific toxicological effects before they are apparent in histopathology studies. In this study, histopathology examination of heart did not show any defect though our method of gene expression analysis could identify enrichment of acute phase response signalling genes in heart that may point towards building up of toxic responses in heart. Given the fact that many PPAR agonist drugs have been shown to cause cardiac toxicity on prolonged usage and FDA's requirement of one year toxicity study for PPAR agonist drugs, our results show promising early detection of toxicity in the drug discovery process.
Overall, our aim here is to establish a proof of principle for the approach of simultaneously analysing multiple, related large datasets. We therefore focused on a dataset where clear-cut validation is possible. However, we note that the technique is very general, and therefore opens up many new opportunities in datadriven computational biology. In particular, it can be applied to heterogeneous sources of data; for example, generated by different laboratories or experimental methodologies. We are currently pursuing this approach in the study of colon cancer.