Investigation of the usefulness of liver-specific deconvolution method toward legacy data utilization

Background Immune responses in the liver are related to the development and progression of liver failure, and precise prediction of their behavior is important. Deconvolution is a methodology for estimating the immune cell proportions from the transcriptome, and it is mainly applied to blood-derived samples and tumor tissues. However, the influence of tissue-specific modeling on the estimation results has rarely been investigated. In this study, we constructed a system to evaluate the performance of the deconvolution method on liver transcriptome data. Results We prepared seven mouse liver injury models using small-molecule compounds with known hepatotoxicity and established a dataset with corresponding liver bulk RNA-Seq and immune cell proportions, covering various immune responses. RNA-Seq expression for nine leukocyte subsets and four liver-associated cell types were obtained from the Gene Expression Omnibus (GEO) to provide a reference covering liver component cells. Here, we found that the combination of reference cell sets affects the estimation results of reference-based deconvolution methods. We established a liver tissue-specific deconvolution by optimizing the reference cell set for each cell to be estimated. We applied this model to independent datasets and showed that liver-specific modeling focusing on reference cell sets is highly extrapolatable. Conclusions We provide an approach of liver-specific modeling when applying reference-based deconvolution to bulk RNA-Seq data and show its importance. It is expected to enable sophisticated estimation from rich tissue data accumulated in public databases and to obtain information on aggregated immune cell trafficking.


Conclusions
We provide an approach of liver-specific modeling when applying reference-based deconvolution to bulk RNA-Seq data and show its importance. It is expected to enable sophisticated estimation from rich tissue data accumulated in public databases and to obtain information on aggregated immune cell trafficking.

Background
The data analysis method called deconvolution extracts immune cell information such as the proportion of immune cells in the sample from the bulk transcriptome data. The bulk transcriptome has a long history, starting with microarrays around 2000, and lots of data have been accumulated in public databases and are easily available [1]. Therefore, it is expected that the combination of deconvolution method with these legacy data enables us to obtain aggregated knowledge on the behavior of various immune cells, so called immune cell trafficking, under various conditions. To obtain such aggregated knowledge, it is difficult to utilize flow cytometry data, which is not yet available as a well-organized database, and single-cell RNA sequencing (scRNA-seq) data, which has been in vogue in recent years but has accumulated relatively little data [2].
Typical deconvolution methods are reference-based methods that use the unique gene expression levels of immune cells as prior information [3][4][5]. When these methods are applied to tissue data accumulated as legacy data, how accurate is the estimation? While evaluation data sets exist for blood, and their performance has been evaluated [3,[6][7][8], is there any tissue dependence when applying these methods to tissues with other parenchymal cells? In fact, such a possibility has been pointed out by Ziyi Chen et al. [9]. However, since no evaluation dataset exists, there was no clear answer.
In the present study, we analyzed mouse liver tissues to verify the accuracy of the reference-based deconvolution method in tissues. Immune cell trafficking, especially neutrophils, is important in the development and progression of liver failure and studies using various liver injury models have reported the role of immune cells in the migration of neutrophils [10][11][12][13][14]. However, the studies have been limited to individual models of specific disorder, and there is no aggregated knowledge of the commonalities and differences in the contribution rate and the order of elicitation of each cell type.
Therefore, we evaluated the accuracy of the deconvolution method by establishing an evaluation dataset covering various immune cell behaviors using compound-induced liver injury models, in which there is relatively little confounding among the models.
The present study has three contributions to the understanding of immune responses in tissues as follows: (1) Using several liver injury models induced by each of seven small-molecule compounds with known hepatotoxicity, we constructed a dataset for evaluation of the deconvolution method covering various immune responses.
(2) In the reference-based deconvolution method, we showed that it is important to select cell types that constitute the reference.
(3) By using a liver-specific optimized model, we were able to find behaviors of immune cells that could not be found by conventional methods.

Establishment of evaluation dataset using drug-induced liver injury models
There is no evaluation dataset that corresponds to the liver sample transcriptome and immune cell type proportions. To construct a liver-specific deconvolution model, we first prepared an evaluation dataset that reflect the immune response in liver tissue. In general, the evaluation dataset should be as diverse as possible in terms of input-output relationships and should be free of confounding factors other than inputs and outputs. We prepared a diverse and small confounding dataset using perturbation by various small molecules for the application to low-throughput in vivo studies. Liver injury models with small molecular compounds are expected to contribute to the establishment of a highly diverse data set, with little confounding by procedures such as surgery. Based on a literature survey, we selected seven compounds, alpha-naphthyl isothiocyanate (ANIT), acetaminophen (APAP), CCl4 (CCl4), concanavalin A (ConA), galactosamine (GAL), 4,4'-methylene dianiline (MDA), and thioacetamide (TAA), that are frequently used to induce liver injury [15][16][17][18][19][20][21].
Liver injury was induced by administration of each compound, and tissues were collected 24 hours after administration. Flow cytometry, RNA-seq analysis, and blood biochemistry tests were performed on the obtained tissues, and evaluation dataset was established in which these measurements corresponded to each sample (Fig. 1a).

RNA-seq analysis on the evaluation dataset
Liver transcriptome data were obtained by RNA-Seq. We performed PCA on the processed data and showed that each treatment group formed a cluster (Fig. 1b). the average expression levels of the genes in each treatment group were calculated and visualized in a heat map, which detected differentially expressed genes (DEG) among the treatment groups (Fig. 1c). These results indicate that each of the seven compounds used to induce liver injury had specific effects on the tissue and that gene expression profiles were separated among the treatment groups. Thus, our established evaluation dataset using druginduced liver injury models is expected to reflect the diverse immune responses in the liver (Fig. S1).

Flow cytometry analysis on the evaluation dataset
The flow cytometry data were analyzed with FlowJo software. In the present study, we obtained ground truth data on the proportion of nine subsets representative in the liver, αβT cells, γδT cells, natural killer T (NKT) cells, natural killer (NK) cells, monocytes, neutrophils, Kupffer cells, eosinophils, and monocytederived macrophages (MDMs), in the compounds induced liver injury samples (Fig. S2). Each immune cell can be evaluated for increase or decrease compared to untreated control samples We observed some characteristic immune cell trafficking for each compound administration (Fig. 1d). For example, APAP administration increased neutrophils, ConA administration increased monocytes, and TAA administration decreased NK cells. These results were consistent with existing findings [22][23][24]. On the other hand, to the best of our knowledge, this is the first report of an increase in eosinophils after GAL administration and the responses of various immune cells after MDA administration captured using flow cytometry.
It is noteworthy that, in this present analysis, each sample is tied to the value of liver injury markers, so we can associate the degree of damage with the behavior of immune cells (Fig. S3). There were large individual differences in liver injury in the APAP and ConA treatment groups, respectively, which could be stratified by the accumulation of immune cells. Specifically, the trafficking of eosinophils after APAP administration varied widely, but there was a tendency for eosinophils not to be accumulated in the severely injured samples with high alanine aminotransferase (ALT) and aspartate aminotransferase (AST).
Furthermore, ALT and AST values at ConA administration correlated with the accumulation of NK cells and Neutrophils ( Fig. 1e and S4). In the present study, the established evaluation dataset was used to evaluate and optimize the deconvolution method, but we would also like to emphasize that it provides novel insights useful for understanding the mechanisms underlying the induction of liver injury and the immune cell trafficking by each compound administration.

Evaluation of the impact of the reference cell sets on estimation performance
The immune cell proportions were estimated by deconvolution using Elastic Net on the established evaluation dataset. In the present study, we obtained cell-specific gene expression profiles for up to 13 types of cells derived from immune cells and liver cells and made them available as reference (supplementary tableS1). The correlation and multicollinearity of gene expression profiles of these 13 cell types are shown in Figure S5.
We prepared three representative reference cell sets. The first is a set of six cell types called LM6 (Neutrophils, Monocytes, B, CD8, CD4, and NK), which are frequently used in estimating the proportion of immune cells in blood using the deconvolution method [25,26]; the second is LM9, to which we added immune cells that are not frequently used but are important in the liver (e.g., Eosinophils, Basophils, and Kupffer cells); and the third is LM13, to which we added four cell types related to the liver. Using these references, we evaluated the estimation performance on the evaluation dataset ( Fig. 2a-c).
For the estimation of neutrophils, LM6 showed superior performance and LM13 showed poor performance. On the contrary, for the estimation of NK cells, LM6 performed poorly, but LM13 performed well. The reference cell types dependence was also confirmed for other representative algorithms other than Elastic Net (Fig. S6). These findings indicate that the estimation accuracy of each cell depends on the combination of cell sets placed in the reference, and that there is no reference with consistently good estimation performance (Fig. 2d). Since the Elastic Net showed the best performance when all 13 cell types were used, we optimized this method for subsequent analyses (Fig. 2e).

Optimization of the combination of reference cell types
In general, deconvolution is a methodology to estimate the composition of cells. On the other hand, when describing individual tissue conditions, such as elucidating disease mechanisms or stratification, it is important to focus on a specific cell and estimate how the cell fluctuates compared to the control group in terms of trafficking. In the case of trafficking, it is possible to build a model for each cell because it is specified for a control group. Therefore, we optimized the model for each target cell and evaluated the characteristics of the reference cell sets.
When estimating the trafficking of a focused cell, all combinations of 12 cell types other than the focused cell are candidates for the cell set to be placed as a reference. For the cell to be analyzed, differentially expressed genes (DEGs) were generated as a reference for all combinations, and the estimation performance of the Elastic Net was evaluated. In this evaluation, the samples in the dataset were those in which the target cells of the trafficking analysis showed fluctuation in comparison to the control group by flow cytometry.
First, the influence of all combinations of cells in the reference was evaluated using the correlation coefficient between the deconvolution estimated value and the actual value measured by flow cytometry for each of neutrophils, monocytes, NK cells, and eosinophils (Fig. 3a). The correlations between the estimated and measured values when the top 10 and bottom 10 combinations were calculated, and differences in estimation in performance were observed depending on the choice of reference cell sets  Fig. 3d and S8). It was shown that the combination of cells to be placed in the reference has good and bad affinity is also existed in the estimation of other cells (Fig. 3e). Therefore, we selected neutrophils, monocytes, NK cells, and eosinophils as the target cells for estimation due to computational cost issues and examined combinations of reference cells to establish a liver-specific optimized deconvolution model ( Fig. S9-S12).

Application to external public data
To aggregate immune cell findings from legacy data based on this model, it is essential to evaluate extrapolation. We then applied our optimized models based on reference cell combinations to external data for a more realistic task and evaluated their extrapolation and robustness. Bird et al. have published transcriptome data at different time points after APAP administration (GSE111828) [27]. Since no immune cell information accompanies this dataset, we applied deconvolution using Elastic Net to each time point and evaluated the differences in performance with and without the liver-specific optimization.
We compared the number of cell types constituting the top 10 and bottom 10 references in the estimation performance of the four cell types investigated for all reference cell combinations: neutrophils, monocytes, NK cells, and eosinophils (Fig. 4a). In preparing the references according to whether they are optimized or not, we attempted to eliminate the confounding effect of the number of cell types placed in the references on estimation performance. The threshold was set as the least number of cell types among the top 10 reference cell combinations found in the optimization process using the evaluation data set. Nonoptimized combinations were selected for the bottom 10 combinations that retained cell types above the threshold. We performed deconvolution on GSE111828 with and without optimization, respectively, and the estimated values were also compared as baseline when all 13 cell types were used. The results showed differences in the estimated percentage of immune cells at each time point after APAP administration ( Fig.   4b). For example, the optimized model estimated an increase in eosinophils with a peak at 48 hours after APAP administration, while the non-optimized model had the opposite estimate of a decrease. The estimates for neutrophils at 48 hours post-dose also showed opposite behaviors of increase and decrease, respectively, with and without optimization.
Next, we evaluated the differences in the GSE111828 immune cell proportion estimates at each time point with and without the optimization described above to see which reflected true immune cell behavior. We performed a reproduction study of APAP administration following the report of Bird et al. and measured immune cell proportions by flow cytometry to evaluate the extrapolation of the optimized model (Fig. 5a).
For lymphocytes, six subtypes were selected as immune cells to be measured using flow cytometry:

CD4+ T cells (CD4), CD8+ T cells (CD8), gamma-delta T cells (γδT), B cells, natural killer T (NKT) cells, and natural killer (NK) cells. For bone marrow-derived cells, measured immune cell proportion data
were obtained for five subsets: monocytes, neutrophils, eosinophils, Kupffer cells, and monocyte-derived macrophages (MonoMac). Gating in flow cytometry of lymphocytes and bone marrow-derived cells is described in Figure S13. Each immune cell proportion at 12, 24, and 48 hours after APAP administration was obtained and showed characteristic changes (Fig. 5b).
The correspondence between the estimated values and the actual values measured by flow cytometry for neutrophils, monocytes, NK cells, and eosinophils at each time point of GSE111828 was evaluated (Fig.   5c). Flow cytometry measurements showed an increase in eosinophils and a decrease in neutrophils at 48 hours after APAP administration. These results were consistent with those estimated using the optimized reference, and the establishment of a tissue-specific deconvolution model enabled precise prediction.
Furthermore, the findings are consistent with recent reports [28]. These results indicate that the combination of references optimized in this study is not merely an overfitting to the prepared evaluation dataset but is highly extrapolative.     The present study also revealed a tendency for eosinophils not to be induced in severely injured samples with high levels of alanine aminotransferase (ALT) after APAP administration. There is no stratification of eosinophils' trafficking according to the degree of liver injury after APAP administration, and it is possible that the molecular mechanisms of hepatic protective function responsible for eosinophils' trafficking may be disrupted during severe injury. These interesting data on eosinophils would contribute to the elucidation of a common molecular mechanism in the hepato-protective function of eosinophils in a variety of liver dysfunctions with different toxic action points and mechanisms.
In the optimization process of the reference cell sets using the established evaluation dataset, there was a correspondence between the biological meaning and the combinations, such as Kupffer cells is frequently included in the combinations suitable for neutrophil estimation. This may be due to the multicollinearity of the deconvolution model. Although the Elastic Net employed in the deconvolution model of this study can handle multicollinearity with L2 penalty to some degree [31], the effect itself is present. CD4+ T cells and CD8+ T cells show similar gene expression profiles and are highly multicollinear (Fig. S5) The present study has several potential limitations. We took the liver as an example of a tissue and showed the optimization of tissue-specific deconvolution and its usefulness. The liver is a single organ with a variety of disorders and was suitable in that it is possible to construct models with less confounding using small molecular compounds [32]. On the other hand, when performing deconvolution for tissues other than the liver, it is necessary to prepare disorder models that reflect various immune cell trafficking to the target tissue and to construct an optimized model using the evaluation dataset. Also, for human clinical tissue samples, it is expected to be difficult to establish evaluation dataset using compound perturbations, and the pipeline in the present study cannot simply be applied. Therefore, it is important to apply deconvolution methods other than reference-based methods that require tissue-specific optimization.

Conclusions
In recent years, tissue transcriptome data have been abundantly accumulated in databases, and the acquisition of novel biological insights from these legacy data is an important and challenging task.
Deconvolution methods that extract immune cell information from tissue transcriptome data are attractive options. Here, for liver-specific modeling, we established an evaluation dataset of the corresponding transcriptome and immune cell proportions in various mouse liver injury models using small compounds.
Using this dataset, we provided a methodology to optimize the combination of reference cell types for each cell to be estimated. The optimized model was applied to external data, suggesting that the model more accurately captures the time-series changes of each immune cell after APAP administration and is highly extrapolatable. The results presented here emphasize the importance of tissue-specific modeling when applying reference-based deconvolution methods to tissue transcriptome data.

Animals
Five-week-old and nine-week-old male C57BL/6JJcl mice were purchased from CLEA Japan, Inc.
(Tokyo, Japan) and kept under standard conditions with a 12-hour day/night cycle and access to food and water ad libitum.

Drug-induced liver injury models
Six-week-old mice acclimatized for one week were used for the experiments. Seven compounds, alpha-  Supplementary Table S2. As a control group, 0.5% methylcellulose solution was used for oral administration, and saline was used for tail vein and intraperitoneal administration. The dose was 10 µL/kg for all groups. Twenty-four hours after administration, animals were sacrificed and perfused liver samples and blood were harvested (Liver and blood sample collection section). For the obtained liver, a portion of the outer left latera lobe of the liver was subjected to RNA isolation (RNA-seq analysis section) and the remaining tissue to flow cytometry analysis (Flow cytometry analysis section), respectively.

Time dependent model of APAP-induced liver injury
Bird et al. performed RNA-seq analysis on mouse liver tissue samples at different timepoints after APAP administration [27], and we followed the protocol.
Ten-week-old mice acclimatized for one week were fasted for 10 hours from 7 am to 5 pm. APAP was dissolved in sterile phosphate-buffered saline (PBS) warmed to 42°C and administered at 350 mg/kg by single i.p. injection of 20 µL/g. Animals were sacrificed 12, 24, and 48 hours after APAP administration, and perfused liver samples and blood were collected (Liver and blood sample collection section). The harvested liver sample was subjected to flow cytometry analysis (Flow cytometry analysis section).

Liver and blood sample collection
Briefly, a superior vena cava was clipped using a clamp and the blood was collected through an inferior vena cava into 1.

Flow cytometry analysis
The liver sample was dissociated using gentleMACS according to the manufacturer's instructions. Except noted especially, phosphate-buffered saline containing 2% fetal bovine serum was used as "Wash Buffer" thereafter. The washed samples were centrifuged (50×g, 4 °C for 3 min) to eliminate hepatocytes and were subjected to ACK Lysis Buffer.

Elastic Net
Consider a measured bulk gene expression matrix ∈ × for N genes across M samples, each containing K different cell types. The goal of deconvolution is to estimate cell type-specific expression ∈ × and cell type proportion matrix ∈ × , and can be written as: Elastic Net [31] is a regularized regression model with combined L1 and L2 penalties. We can estimate the cell type proportion matrix via: where and ! are hyper parameters, and we set = 1 and ! = 0.05 as default parameters.

Bulk tissue gene expression matrix
TTPM normalized liver expression profile is the analysis target. The transcript IDs were converted to MGI gene symbols using files available from Biomart [34], and median values were selected for duplicate gene names. The data were processed in the following order: log transformation, elimination of lowexpressing genes, batch normalization, and quantile normalization. Among the genes in this normalized expression matrix, we focused on N marker genes, which are described below, and defined them as ∈ × .

Cell type-specific expression matrix (Reference)
First, we downloaded raw RNA-Seq expression datasets for 9 leukocyte subsets and 4 liver related cell types from the Gene Expression Omnibus (GEO) Hierarchical clustering with Pearson correlation coefficient was performed on TPM normalized profiles, and samples forming the main cluster were manually selected. Note that data is available at our github repository (refer to Availability of data and materials section). Then, we converted transcript IDs to MGI gene symbols and median values were selected for duplicate gene names. The data were processed in the following order: log transformation, elimination of low-expressing genes, and quantile normalization. Only genes in common with the bulk tissue gene expression matrix were retained. Among the retained ones, up to 50 genes were kept as markers with absolute fold change greater than 1.5 for the second cell type with highest expression.
Finally, N genes collected from K cell markers defined above were included in the analysis, and cell typespecific expression matrix ∈ × was defined.

Evaluation
Estimated cell type proportion matrix ∈ × was obtained by running deconvolution and converted to a sample-wide z-score for each cell of interest. Similarly, the ground truth cell type proportion matrix obtained by flow cytometry was converted to a sample-wide z-score and its Pearson Correlation with was evaluated. Note that the estimated is not necessarily greater than 0 because it is not constrained to be nonnegative.

Selection of samples showing fluctuation
We

Ethics approval and consent to participate
The studies reported in this article were performed in accordance with the guidelines provided by the Institutional Animal Care Committee (Graduate School of Pharmaceutical Sciences, the University of Tokyo, Tokyo, Japan).

Consent for publication
Not applicable.

Availability of data and materials
Code, models, and data are available at https://github.com/mizuno-group/LiverDeconv. RNA-seq data will be shown in GEO datasets.

Competing interests
The authors declare that they have no conflicts of interest.             Table S1. Summary of sample information used to construct a reference consisting of up to 13 cell types.