Detecting networks of genes associated with human drug induced liver injury (DILI) concern using sparse principal components

Abstract

The 12th Annual International Conference on the Critical Assessment of Massive Data Analysis (CAMDA) used data from the massive Japanese Toxicogenomics Project (TGP) to predict drug-induced liver injury (DILI) concern provided by the U.S. Food and Drug Administration (FDA). The challenge was to predict DILI concern by means of gene expression data. Analysis of this high-dimensional toxicogenomic data requires statistical methodologies that can detect the transcriptomic associations with toxicity. We propose an analysis technique that involves sparse principal component analysis to efficiently reduce the dimension of the analysis problem. Sparse principal component variables are composed of groups of expressed genes. Associations between DILI concern and sparse principal component variables were tested and further scrutinized with sparse regression methodology to identify concise transcriptomic structures potentially responsible for and predictive of drug toxicity. Working with a subset of the TGP data with FDA DILI concern classification, we identified 5 transcriptomic structures (sparse principal component variables) statistically associated with DILI concern. The most statistically significant structure consists of the genes ZBTB16, FLVCR2, TNS3, and ASB13. Sparse statistical methods offer a new way to handle analysis issues with massive omic data. Sparse PCA can efficiently extract groups of transcriptomic markers that may indicate drug toxicity.

Keywords:

• statistics
• toxicogenomics
• drug toxicity
• drug-induced liver injury
• sparse principal component analysis
• sparse principal components
• loadings
• groups of genes
• differentially expressed

Abbreviations

CAMDA	=	Critical Assessment of Massive Data Analysis
DEG	=	differentially expressed gene
DEPC	=	differentially expressed principal component
DILI	=	drug-induced liver injury
FARMS	=	Factor Analysis for Robust Microarray Summarization
PC	=	principal component
PCA	=	principal component analysis
TGP	=	Japanese Toxicogenomics Project

Introduction

The current estimated cost in US dollars to develop a drug is $1.8 billion, but most drugs never make it to market, largely due to toxicity levels deemed unsafe for the human liver.¹ As it is unethical to test the effects of compounds in humans, toxicity has classically been tested in animal-based experiments. Hence, conclusions may not be generalizable to the human population, especially when the effects of toxicity are revealed after prolonged exposure. The testing inaccuracy could result in passing unsafe or discarding useful drugs, which hinder public health efforts and squander resources. In an attempt to improve testing accuracy, the field of toxicology has embraced animal and now human “omic” data (genomics, transcriptomics, proteomics, metabolomics) to help identify biological material that predict the toxicity of compounds at an early stage. Therefore, resources are saved from developing drugs that will ultimately fail in the public domain. The use of omic data to inform toxicology is commonly referred to as toxicogenomics and is now the focus of several research initiatives. One such initiative is the Japanese Toxicogenomics Project (TGP).² With human in vitro, rat in vitro, and rat in vivo experimental models, 131 compounds were applied to liver samples and microarray gene expression profiles were obtained using Affymetrix GeneChip® technology.

To facilitate prediction of toxicity in new drugs, the Food and Drug Administration (FDA) developed a classification system of human drug-induced liver injury (DILI) concern (most, less, and no DILI concern) for drugs currently on the market.³ The drugs they classified had been on market for a minimum of 10 y, allowing sufficient public interaction to obtain updated and realistic DILI concern information, not otherwise attainable due to ethical reasons. This new classification of human-based drug toxicity can be linked to toxicogenomic databases, facilitating the search for omic markers that predict toxicity based on this indicator. The 12th Annual International Conference on Critical Assessment of Massive Data Analysis (CAMDA 2013) linked FDA classification labeling to the TGP data and proposed analysis challenges involving prediction of drug toxicity. Discovering novel biomarkers associated with DILI concern within the TGP data may aid in our understanding of mechanisms of toxicity and enhance our ability to assess the toxicity of new compounds. However, the breadth and complexity of omic data makes analysis challenging and to extract key information it is important to enrich statistical methodology with biological context.⁴

Commonly, simple statistical methods are used to test associations between drug toxicity and each marker, one at a time, developing a list of top candidate genes. Although this approach is easy to implement, it involves conducting thousands of statistical tests and is prone to spurious associations (i.e., the “multiple testing problem”). Perhaps more concerning is that by testing markers independent of one-another, this analysis does not acknowledge the fact that genes may act in concert to influence a phenotype. Complex networks of omic material will likely be the underlying indicator of drug toxicity. Identifying and characterizing these indicators requires more sophisticated statistical methods. We believe that recent advancement in a new class of so-called sparse statistical methods validates the use of principal component analysis (PCA) as a primary analysis tool to detect these complex networks.

PCA is a commonly used multivariate method for both dimension reduction and data visualization. It assembles a new set of variables, called principal components (PCs), from linear combinations of original variables. These PCs are uncorrelated and ordered by maximal variance, giving the analyst an easier data set to work with if they can interpret what the PCs mean. Unfortunately, this is a major challenge with high-dimensional data, as found in toxicogenomics, since the PCs are formed by linear combinations of all original variables (a weighted sum of 1000 genes, for example) and are uninterpretable. To overcome this limitation, sparse principal component analysis methods^5-7 have been developed by combining PCA with sparse regression methodology.^8,9 Sparse PCA restricts principal components to be formed by interpretable linear combinations of smaller (sparse) subsets of the original variables (a weighted sum of 10 genes, for example). Tuning parameters control the level of sparseness induced, making sparse PCA procedures very flexible. Sparse principal components reflect concise and interpretable groups of original variables that remain in the linear combinations.

In this paper, we present an analysis strategy built around using sparse PCA to extract groups of related genes from a toxicogenomic database, testing their associations with drug toxicity. We apply this analysis strategy to subsets of the TGP database with the FDA's classification of human DILI concern as the measure of drug toxicity. We hope to convey that this new class of sparse statistical methods may prove beneficial to toxicogenomic analyses, as they specialize at separating the “signal” from “noise” amidst high-dimensional data.

The rest of our paper is organized as follows. In the Materials and Methods section, we describe the data we used, the essentials of sparse methodology, and our analysis strategy. We then present our Results and provide a Discussion on the potential merits and limitations of sparse methodology in toxicogenomic analyses.

Materials and Methods

Description of Data

Data source

Completed in 2007, the TGP database was the result of a 5-y collaborative effort between government and private companies to obtain gene expression data using Affymetrix GeneChip® technology to characterize the response to various drugs across a variety of experimental units.² The completed TGP database provided by CAMDA 2013 contains thousands of results from experiments, involving 131 drugs applied at 4 dose levels (control, low, middle, and high) to human and rat in vitro hepatocytes, and administered to rat in vivo with either a single dose or repeated doses. Transcriptomic material was extracted from experimental units at several time points after the drugs had been given. Full details regarding the study design and protocol are given by Uehara et al.² and complete data, along with the portion we obtained including FDA human DILI concern labels, are available through the CAMDA 2013 conference website (http://dokuwiki.bioinf.jku.at/doku.php).

Samples used

We considered only those drugs with FDA human DILI concern classification; 93 of the 119 drugs tested on human samples, and 101 of the 131 drugs tested on rat samples. We primarily focused on the 93 samples from human in vitro experiments that specifically received high dose levels and had gene expression measured at 8 h; a single subset of this rich database. However, to compare results across different experimental conditions, we also repeated our analysis on another 15 subsets of the database. Anticipating that higher doses result in more robust gene expression measurements¹⁰ and due to many drugs not being administered at low doses in the human in vitro samples, we considered only samples that received middle or high dose levels. Likewise, since measurements of gene expression in the human samples were not taken at 2 h for many drugs, we considered only gene expression values measured at later time points.

Transcriptomic markers used

Human gene expression data were obtained with the Affymetrix GeneChip® Human Genome U133 Plus 2.0 Array and rat gene expression data were obtained with Affymetrix GeneChip® RAE 230A 2.0 Array.² Although raw gene expression data was available, we chose to use preprocessed data that had undergone batch-effect correction and the Factor Analysis for Robust Microarray Summarization (FARMS),¹¹ that was also provided. A total of 18988 probesets were available for replicate-collapsed human samples, and 12088 probesets were available for replicate-collapsed rat samples. Almost all probesets had gene names provided and all gene names were unique; as such we refer to probesets and genes interchangeably unless the gene name was not provided. We used inter-quartile range to filter genes that did not vary much and kept the top 1000 genes to simplify analysis.

Classification variable for human DILI concern

Human DILI concern (“most,” “less,” and “no” DILI concern in humans) was provided to detect if gene expression measurements differed across DILI classification. However, since only 8 drugs are classified as “no DILI concern,” we reclassified the human DILI concern variable to be binary (“most” vs. “less or no” DILI concern) to obtain a more balanced and simple categorical variable. Of the 93 drugs tested in human tissue, 40 are labeled as “most DILI concern” and 53 are “less or no DII concern” and of the 101 drugs tested in rat specimens, 41 are “most DILI concern” and 60 are “less or no DILI concern.”

Sparse PCA and Sparse Regression Methodology

Here we give an intuitive description of sparse PCA and sparse regression and refer the reader to refs.^5-9 for mathematical details.

Given a data matrix X = (X₁,X₂, …, X_p), where p variables (e.g., probsets or other omic markers) are measured for n observations (samples; drugs), classical PCA constructs linear combinations (weighted sums) from the variables in the form${\displaystyle Z_i=v_{i,1}{X}_1+v_{i,2}{X}_2+v_{i,3}{X}_3+\dots +v_{i,p-1}{X}_{p-1}+v_{i,p}{X}_{p,}i=1,2,\dots ,\min \left(n,p\right)}$

where Z_i is the i^th principal component (PC) variable and the v_i,j's are the loadings (coefficients; weights) for the j^th variable in the i^th PC. By utilizing the correlation structure among X's, PCA constructs PCs in such a way that they are uncorrelated, ordered by maximal variance, and the total amount of information (variance) in the new data matrix Z = (Z₁,Z₂, …, Z_min(n,p)) is equal to the total amount of information in the original data matrix X. We may now choose to analyze the Z's instead of X's, as these defining properties potentially simplify analysis. For example, analyzing uncorrelated Z's allows for easier regression analysis as there is no issues with colinearity. Even more appealing, if n < p then the number of Z's to analyze will be less than the number of X's. Likewise, since Z's are ordered by variance, the first few PCs will typically hold a large portion of the information and we can discard the trailing ones without losing much. For these last two reasons, PCA is an excellent so-called dimension-reduction tool, able to compress large amounts of data into a few important components. In some cases, PCA can be an excellent exploratory tool, since groups of highly correlated X's will bear substantial weight in the same PC and we can detect groups of variables through the magnitude of their loadings (v_i,j's). However, since PCs are linear combinations of all X's, this makes them nearly impossible to interpret when dealing with thousands of correlated X's, as is common with high-dimensional toxicogenomics data.

Sparse PCA is an extension to classical PCA that forces the loadings for some X's to be exactly 0, creating sparse linear combinations in the form (for example):${\displaystyle Z_i=v_{i,1}{X}_1+v_{i,2}{X}_2+0X_3+\dots +0X_{p-1}+v_{i,p}{X}_p,i=1,2,\dots ,\min \left(n,p\right)}$

where now the v_i,j's that do not have much relation to a PC are removed. This results in a more interpretable set of PCs, having small groups of X's (genes, for example) creating new variables to work with. This is objectively achieved through integrating constrained (sparse) regression methodology to the derivation of PCA. For example, finding classical PCA solutions through the singular value decomposition (SVD) X = UDV′, where V holds all the loadings and Z = UD are the PCs, Witten et al.⁶ applied constraints to elements of V, shrinking some of them to 0. A so-called tuning parameter (λ) controls the number of loadings that are forced to 0, providing a flexible framework from which to obtain sparse PCs. As a trade-off for acquiring interpretable PCs by introducing sparseness, a proportion of the total information (variance) is lost.

Finally, sparse regression works in a similar way, in that the coefficients estimates (β's) of a regression model of the form${\displaystyle Y={\hat{\beta }}_0+{\hat{\beta }}_1{X}_1+{\hat{\beta }}_2{X}_2+{\hat{\beta }}_3{X}_3+\dots +{\hat{\beta }}_{p-1}{X}_{p-1}+{\hat{\beta }}_p{X}_p+\epsilon ,}$

are constrained, shrinking a number of them, dictated by tuning parameters, directly to 0. It has been integrated to many extensions of the linear model, such as generalized linear models, including binary logistic regression.

Analysis Strategy

The novel component of our analysis strategy is the use of sparse PCA to automatically select groups of variables, however we suggest a more broad analysis pipeline as visualized in Figure 1.

Figure 1. A visual of our analysis strategy applied to human in vitro, high dose, 8 h sampling time data. We begin in the top row (1-a), by conducting a Differentially Expressed Genes (DEGs) analysis on the gene expression matrix; columns represent the 93 samples (40 “most” and 53 “less or no” DILI concern), rows represent 1000 expressed genes. This returns a list of top DEGs (1-b); the genes that are most significantly associated with DILI concern. We then move to the middle row (2-a), using sparse PCA on the same gene expression matrix to obtain new sparse principal component (PC) variables (2-b) to work with; columns for this new data matrix again represent the 93 samples, but rows represent the 93 new sparse PC variables (we have reduced the dimension from 1000 to 93). Then, we conduct a Differentially Expressed PCs (DEPCs) analysis on the PC expression matrix to obtain a list of top DEPCs (2-c); the sparse PCs that are most significantly associated with DILI concern. At this point, we examine the genes that contribute to these DEPCs to makes sense of what the structures mean and make note of those genes in these structures that were also identified as differentially expressed in the DEGs analysis. As a final validation step (3-a), we apply sparse regression to the same 93 sparse PC variables to identify a concise list of sparse PCs that are potentially related to DILI concern (3-b).

We implemented this analysis for each of the 16 subsets of the TGP data we considered but refer only to the human in vitro, high dose level, 8 h gene expression sampling time subset while explaining the steps in detail. All analysis was conducted with R statistical software (version 3.0.1).¹²

Step 1: Identify genes independently associated with DILI concern with a Differentially Expressed Gene (DEG) analysis

This step can be visualized with 1-a and 1-b in Figure 1. Starting with the 1000 most variable genes, we first investigated if any genes are independently associated with DILI concern. Instead of using the two independent samples t test for each gene, we used the moderated t test¹³ from the R package “limma”¹⁴ which yields more conservative P; without this correction, when conducting many tests, suspiciously small standard errors can inflate test statistics. We obtained “moderated t” test statistics with associated P and claimed that genes with P < 0.05 were significantly associated with DILI concern. Additionally, using “limma,” we obtained q values (P adjusted to account for false-discovery rates)^15,16 to eliminate genes likely to be claimed falsely significant due to conducting multiple tests. With this list of differentially expressed genes likely to be independently associated with DILI concern, we moved to search for more complex relationships between the gene expression data and DILI concern.

Step 2: Identify groups of genes jointly associated with DILI concern using sparse PCA and a Differentially Expressed Principal Component (DEPC) analysis

This step can be visualized with 2-a, 2-b, and 2-c in Figure 1. Starting with the 1000 most variable genes, we used the sparse PCA method by Witten et al.⁶ executed with the R package “PMA”¹⁷, to obtain sparse principal component variables. We chose this sparse PCA method over certain other formulations^5,7 as a result from our previous work¹⁸ that found it to perform better in computer simulations with high-dimensional data. For each of a range of tuning parameters (λ = 2, 3, 5, 7, 10, 15), we ran sparse PCA and examined the trade-off between adjusted percentage explained variance¹⁸ and sparseness among PCs with the goal to select a tuning parameter that delivered principal components with extremely sparse loading vectors while keeping a large proportion of variance. We chose tuning parameter λ = 3 because while compared with λ = 5, 7, 10, and 15 it resulted in loading vectors that were drastically more sparse at a small cost of information loss and loading vectors were almost as sparse as those from λ = 2 while comparatively keeping a substantial amount of information. Additionally, and important from a practical view, this delivered small and interpretable linear combinations so we could look at the principal components in clear detail. We then statistically tested associations between each sparse PC variable and DILI concern using a permutation approach claiming sparse PCs with P < 0.05 to be significantly associated with DILI concern. With this list of “differentially expressed” sparse PCs (DEPCs) likely to be independently associated with DILI concern, we then recorded and examined the genes that made up the PCs along with their respective loadings in the PC linear combination.

Step 3: Validate DEPCs with Sparse Regression

This step can be visualized with 3-a and 3-b in Figure 1. We simultaneously entered the sparse principal component variables into a binary logistic regression model framework with DILI concern as the outcome, where the coefficient estimates were obtained with sparse regression methodology; penalized maximum likelihood with a least absolute shrinkage and selection operator (LASSO) penalty⁸ with tuning parameter selected via a cross validation procedure that maintains good prediction from the model, executed with the R package “glmnet.”¹⁹ From this sparse regression model, we obtained the list of sparse PC variables that remained, expecting overlap with the list of sparse PCs identified to be significantly related to DILI concern in Step 2, then recorded and examined the genes and loadings in these PCs. Sparse PCs found both in Step 2 and Step 3 were deemed to be most likely to be associated with DILI concern.

Results

We now present the results from our analysis strategy, reporting details on the human in vitro, high dose, 8 h gene expression sampling time subset but also including summaries of findings for the other 15 subsets.

Of the 1000 genes, 54 had significantly different expression between “most” and “less or no” DILI concern samples, according to the moderated t tests with P < 0.05. Table 1 displays their gene names, effect size, and associated P values.

Table 1. Differentially expressed genes (DEGs; independently associated with DILI concern) from our analysis of the human in vitro, high dose, 8 h gene expression sampling time subset

Rank	Gene name	Diff. in means^a	P^b	Rank	Gene name	Diff. in means^a	P^b
1	SNAPC1	0.228	0.0001	28	TNFRSF1B	−0.096	0.0263
2	TSLP	0.221	0.0005	29	BBS12	0.109	0.0268
3	ANKRD1	0.277	0.0027	30	FKBP5	−0.155	0.0271
4	FHL2	0.191	0.0050	31	MIR22HG	−0.189	0.0273
5	CTH	0.126	0.0055	32	FSTL1	−0.156	0.0280
6	HEXIM1	−0.178	0.0057	33	POR	−0.088	0.0282
7	NFIL3	0.195	0.0076	34	TUFT1	0.153	0.0289
8	ETS1	0.117	0.0078	35	LPIN2	−0.098	0.0304
9	MIR3682	−0.130	0.0090	36	RBMXL1	−0.092	0.0306
10	TUBE1	0.157	0.0098	37	CEBPD	−0.124	0.0330
11	BLNK	0.147	0.0117	38	EXT1	0.106	0.0331
12	CIDEC	−0.095	0.0132	39	SLC7A5	0.122	0.0336
13	HAS3	−0.092	0.0134	40	LMCD1	0.187	0.0337
14	ANGPTL4	−0.401	0.0163	41	ERBB3	−0.107	0.0352
15	ELL2	−0.093	0.0173	42	PDLIM5	0.066	0.0352
16	SRSF6	−0.104	0.0179	43	PDK4	−0.138	0.0352
17	INHBE	0.197	0.0190	44	MT1F	0.072	0.0361
18	TXNIP	−0.137	0.0202	45	GEM	0.206	0.0378
19	F3	0.240	0.0205	46	RBMX	−0.111	0.0380
20	FOXA1	−0.134	0.0209	47	ZBTB43	0.127	0.0381
21	MTHFD2	0.173	0.0218	48	HHEX	−0.139	0.0385
22	SLC25A20	−0.125	0.0220	49	TMEM158	0.074	0.0439
23	HSD17B2	−0.136	0.0240	50	ASNS	0.110	0.0446
24	FOXQ1	−0.134	0.0242	51	ATP8B1	0.095	0.0458
25	DUSP6	−0.125	0.0244	52	TNFAIP3	0.147	0.0464
26	CHAC1	0.125	0.0258	53	C11orf96	0.124	0.0474
27	SERTAD2	0.141	0.0261	54	RTP3	0.106	0.0496

^aDifference in Means is “most” – “less or no” DILI concern; a positive value indicates a larger gene expression value for the most DILI concern group. ^bP indicates the P value obtained from the moderated t test.

None of these 54 genes remained statistically significant when using the FDR-adjusted P value (q value) with cut-off q < 0.10. Table 2 displays counts of the number of DEGs found in each of the 16 subsets we analyzed.

Table 2. Counts of DEGs across all 16 subsets analyzed

Data source	Dose level	Gene expression sampling time	# Genes P < 0.05 (mod. t) total (upreg) [downreg]	# Genes q < 0.10 (FDR) total
Human in vitro	High	8 h	54 (29) [25]	0
Human in vitro	High	24 h	26 (12) [14]	0
Human in vitro	Middle	8 h	74 (32) [42]	0
Human in vitro	Middle	24 h	2 (2) [0]	0
Rat in vitro	High	8 h	29 (8) [21]	0
Rat in vitro	High	24 h	153 (1) [152]	17
Rat in vitro	Middle	8 h	42 (39) [3]	0
Rat in vitro	Middle	24 h	34 (20) [14]	0
Rat in vitro (s.Dose^a)	High	9 h	53 (36) [17]	0
Rat in vivo (s.Dose)	High	24 h	82 (49) [33]	0
Rat in vivo (s.Dose)	Middle	9 h	228 (193) [35]	202
Rat in vivo (s.Dose)	Middle	24 h	140 (133) [7]	5
Rat in vivo (r.Dose^b)	High	15 d	50 (27) [23]	0
Rat in vivo (r.Dose)	High	29 d	108 (21) [87]	1
Rat in vivo (r.Dose)	Middle	15 d	53 (31) [22]	0
Rat in vivo (r.Dose)	Middle	29 d	103 (69) [34]	0

The total number of DEGs, along with the number of those that are upregulated (“most DILI concern” has larger gene expression than “less or no DILI concern”) in brackets and those that are downregulated in square brackets. ^a”s.Dose” means rats received only a single dose of drug. ^b”r.Dose” means rats received repeated doses over time.

The rat in vivo samples receiving single dose at middle dose levels with gene expression measured at 9 h had a notably large amount of DEGs compared with the rest; 228 with 202 remaining statistically significant after the more strict FDR-adjustment.

Applying sparse PCA with tuning parameter λ = 3 to the gene expression data, we obtained the 93 [min(n = 93, p = 1000)] sparse principal components. The minimum, median, and maximum number of genes that contributed to a sparse principal component with non-zero loadings was 11, 18, and 33, respectively. Compared with classical PCA, which would return principal components built from all 1000 genes, this is a substantial amount of sparseness and, therefore, interpretability gained. As a sacrifice for simplifying the data in this way, we lost 29.97% of total variance (information) contained in the original 1000 genes. Considering that sparse PCA simplified our focus to a median of just 18 (1.8% of the 1000) genes per component, therefore giving us groups of genes to explore, and reduced our analysis burden from 1000 gene tests to 93 principal component tests, we consider this a massive gain for a comparably small loss. Similar data reduction statistics were observed for the remaining 15 subsets.

Of the 93 sparse principal component variables, 5 had significantly different cumulative expression values between “most” and “less or no” DILI concern samples, according to the permutation test with P < 0.05. Table 3 provides a summary of these DEPCs, including P, a list of contributing genes, and how many of the contributing genes were found in the DEG analysis.

Table 3. Differentially expressed PCs (DEPCs; associated with DILI concern) from our analysis of the human in vitro, high dose, 8 h gene expression sampling time subset

PC	P value	# of contributing probesets (# that are DEGs)	Contributing genes in order of absolute loading value
PC49	0.0074	12 (1)	ZBTB16, FLVCR2, TNS3, ASB13, MLPH, CUX2, FKBP5, SHROOM3, PANK1, TGFBR2, MYLK, ORC6
PC15	0.0213	16 (0)	PCK1, ID1, RRAGC, HSPA1B, SOWAHC, NPC1, PLEKHB2, RNF19B, STX3, CCNE2, KANSL1, PIM3, FNIP2, GNA13, OTUD1, IRS2
PC13	0.0310	13 (8)	DDIT4, MTHFD2, FGF21, DDIT3, INHBE, TRIB3, TUBE1, HSPA13, CHAC1, SLC7A5, TSLP, ASNS, CTH
PC7	0.0382	15 (2)	NUAK2, F3, C8orf4, ENC1, EDN1, CCL2, LMCD1, PLK2, G0S2, KRT7, C1orf63, TRIM6, FILIP1L, THBS1, FASTKD3
PC72	0.0389	18 (3)	SLC40A1, SNAI2, KRCC1, IRS2, LPIN2, CYP1A1, CBLB, BCL6, HMGB2, CEBPD, ERBB3, HECA, RB1CC1, GNA13, EFNA1, TBC1D8, C6orf203

Also provided are counts of probesets that contribute to respective PCs, along with the number that were individually associated with DILI concern (DEGs) and gene names. P values were obtained from permutation tests.

Figure 2 provides a visual representation of the top 3 of these PCs, showing the different types of groups of genes that are associated with DILI concern, in terms of size and composition.

Figure 2. A visual display for the top 3 differentially expressed (most associated with DILI concern) sparse principal components (DEPCs) from our analysis of the human in vitro, high dose, 8 h gene expression sampling time subset. Larger central circles represent the principal components. Attached to each are the genes that form the linear combinations; probesets (gene names) and loading values are inside the outer circles. Shaded circles represent genes that were found to be independently associated with DILI concern (DEGs), whereas non-shaded circles contain genes that were not. PC15 might bring forth a network of trancriptomic material that is associated with DILI concern, not otherwise being found with more simple statistical tests. PC13 shows us that some marginally associated genes behave similarly.

PC49 consists predominantly of probesets (genes), 7704_at (ZBTB16), 55640_at (FLVCR2), 64759_at (TNS3), and 79754_at (ASB13), with lesser contribution from eight more correlated genes, meaning these may have a joint relationship with drug toxicity. One of these eight, 2289_at (FKBP5), was identified as significantly associated with DILI concern in the DEG analysis. PC15 is a slightly larger network comprised entirely of genes which were not identified in the DEG analysis, meaning these genes might only be identified when accumulated into a composite variable like this sparse principal component. PC15 is driven by probesets (genes) 5105_at (PCK1), 3397_at (ID1), and 64121_at (RRAGC). In contrast, 8 out of 13 genes influencing PC13 had independent associations with DILI concern. Its top 4 weighted probesets (genes), 54541_at (DDIT4), 10797_at (MTHFD2), 26291_at (FGF21), and 1649_at (DDIT3), had only one DEG, meaning these genes are likely related to drug toxicity although they did not independently surface. Table 4 reports the number of sparse PCs that were significantly related to DILI concern across all 16 subsets of data we analyzed.

Table 4. Counts of DEPCs across all 16 subsets analyzed. The total number of DEPCs, along with the number of those which are upregulated (“most DILI concern” has larger PC cumulative expression values than “Less or No DILI concern”) in brackets and those which are downregulated in square brackets

Data source	Dose level	Gene expression sampling time	# PCs P < 0.05 (perm.) total (upreg) [downreg]
Human in vitro	High	8 h	5 (3) [2]
Human in vitro	High	24 h	1 (1) [0]
Human in vitro	Middle	8 h	12 (8) [4]
Human in vitro	Middle	24 h	0 (0) [0]
Rat in vitro	High	8 h	2 (2) [0]
Rat in vitro	High	24 h	18 (8) [10]
Rat in vitro	Middle	8 h	3 (1) [2]
Rat in vitro	Middle	24 h	2 (1) [1]
Rat in vivo (s.Dose^a)	High	9 h	3 (1) [2]
Rat in vivo (s.Dose)	High	24 h	5 (3) [2]
Rat in vivo (s.Dose)	Middle	9 h	15 (5) [10]
Rat in vivo (s.Dose)	Middle	24 h	19 (12) [7]
Rat in vivo (r.Dose^b)	High	15 d	6 (4) [2]
Rat in vivo (r.Dose)	High	29 d	12 (6) [6]
Rat in vivo (r.Dose)	Middle	15 d	6 (3) [3]
Rat in vivo (r.Dose)	Middle	29 d	12 (5) [7]

^a”s.Dose” means rats received only a single dose of drug. ^b”r.Dose” means rats received repeated doses over time.

Lastly, sparse regression was used to narrow down a more concise list of DEPCs. Entering the 93 sparse PCs, the model selected PC49 and PC15. These are the top 2 sparse PCs we identified to be most significantly associated with human DILI concern in the DEPCs analysis, bringing some statistical validation to these findings.

Discussion

In this paper, we identified sparse PCs that are differentially expressed between groups of “most” and “less or no” DILI concern, providing small networks of genes for further study in relation to drug toxicity. Apart from being able to investigate groups of variables, there are several components that make sparse PCA attractive. Since the number of principal components returned is min(n,p) (i.e., minimum of n samples and p variables), it will always reduce the number of variables to analyze when applied to high-dimensional (n < p) data. This is greatly beneficial for n < < p data, as we observed with the TGP data set; reducing from 1000 to 93 tests. By reducing the dimensional complexity of the analysis this strategy also reduces the burden of multiple testing, but it does not eradicate the issue and necessary adjustments should be made before concluding if PCs are statistically significant. In contrast to classical PCs, we argue that sparse PCs are interpretable since they are linear combinations consisting of small subsets of genes; however biological context is still required to interpret their values. Built from gene expression data, the sparse PCs in this paper are weighted sums of gene expression values, meaning they can be interpreted as accumulated gene expression across the genes involved. Within the same PC, genes with loading values that have opposite signs are negatively correlated, giving additional insight into structures between genes in the same PC. Since our method returns composite variables, as opposed to simply reporting groups of similar genes as with cluster analysis,²⁰ this allows for statistical testing of joint associations. Since PCA methods are unsupervised (meaning they are not informed by the outcome of interest, such as DILI concern), it is reasonable to expect to find PCs built exclusively from genes that are not differentially expressed. Any of these types of group structures that are found to be associated with the variable of interest potentially offer a brand new source of transcriptomic material for researchers to explore. PC15 in our analysis of the TGP data is an example of this.

We have some cautionary notes regarding the selection of tuning parameters when using sparse PCA to simplify data in search of groups of genes, as it can have major effects on characteristics of the sparse principal components obtained. By choosing λ = 3, for example, we committed to identifying groups of genes of a specific size; a median of 18 genes contributed to building the sparse principal components we obtained. For the purpose of demonstrating immense sparseness, this tuning parameter was a good choice. However, perhaps the true transcriptomic structure underlying the TGP gene expression data involves hundreds of correlated genes expressing similarly and as such, perhaps choosing a less sparse tuning parameter would more accurately capture this truth. This subjective choice of tuning parameter is analogous to restricting the number of clusters, k, in the k-means clustering algorithm. Prediction-based cross-validation methods are available to guide the choice of tuning parameters⁶ and are typically the default methods built into packages.¹⁷ Bonner¹⁸ suggested to visually inspect the variance-covariance structure of data to observe blocking of variables (genes), but this strategy has yet to be tested with rigour or under realistic data conditions.

The analysis strategy we proposed is a practical way to explore multivariate patterns in large genomic data and suggest if they are worth pursuing. Sparse methodology has many extensions to a variety of data analysis situations and it has potential to benefit research in toxicogenomics. The R scripts and analysis pipeline we used in this paper are available upon request.

Disclosure of Potential Conflicts of Interest

No potential conflicts of interest were disclosed.

Funding

JB would like to acknowledge funding from the Natural Sciences and Engineering Research Council of Canada (NSERC) and Canadian Institutes of Health Research (CIHR) (grant number 84392). JB is the inaugural holder of the John D. Cameron Endowed Chair in the Genetic Determinants of Chronic Diseases, Department of Clinical Epidemiology and Biostatistics, McMaster University.