Gene Expression Predicts Histological Severity and Reveals Distinct Molecular Profiles of Nonalcoholic Fatty Liver Disease

Abstract

The heterogeneity of biological processes driving the severity of nonalcoholic fatty liver disease (NAFLD) as reflected in the transcriptome and the relationship between the pathways involved are not well established. Well-defined associations between gene expression profiles and disease progression would benefit efforts to develop novel therapies and to understand disease heterogeneity. We analyzed hepatic gene expression in controls and a cohort with the full histological spectrum of NAFLD. Protein-protein interaction and gene set variation analysis revealed distinct sets of coordinately regulated genes and pathways whose expression progressively change over the course of the disease. The progressive nature of these changes enabled us to develop a framework for calculating a disease progression score for individual genes. We show that, in aggregate, these scores correlate strongly with histological measures of disease progression and can thus themselves serve as a proxy for severity. Furthermore, we demonstrate that the expression levels of a small number of genes (~20) can be used to infer disease severity. Finally, we show that patient subgroups can be distinguished by the relative distribution of gene-level scores in specific gene sets. While future work is required to identify the specific disease characteristics that correspond to patient clusters identified on this basis, this work provides a general framework for the use of high-content molecular profiling to identify NAFLD patient subgroups.

Figure 1

Integration of differentially expressed genes with a protein-protein interaction network highlights hubs involved in the progression of fatty liver disease. (A,B) A protein-protein interaction network induced by the differentially expressed genes for both NAS and fibrosis stage, respectively. Each node represents a densely connected community of proteins, whose size represents the number of proteins in the community. The node labels provide a summary of the biological processes enriched in each community, as well as a number which is a community identifier. Edge thickness is proportional to the number of connections between communities. (C,D) Box plots showing the distribution of eigenvalue centrality in the communities of each network. Communities significantly enriched with hubs (nodes with relatively large centrality) are labeled with their top 5 genes by centrality.

Figure 2

The top Reactome gene sets that are up- and down-regulated with respect to NAS (A) or fibrosis stage (B). The y-axes represent the GSVA score, which is a pathway-level quantification of gene abundance, and the x-axes represent the clinical assessment. For disease activity (NAS), pathways related to apoptosis, inflammation (Fc epsilon receptor signaling, TNFR2 signaling, T cell receptor (TCR)), cell proliferation (PTEN, TP53) were top pathways whereas for insulin receptor substrate (IRS) signaling pathway was downregulated. For fibrosis, Ephrin signaling related genes were the top pathway while amine derived hormones and nicotinic acetylcholine receptor pathways were down-regulated.

Figure 3

Based on the dynamic range of expression and rank order upon ordinal regression of gene expression levels to the NAFLD activity score (NAS) or fibrosis stage, a gene-level score was derived for all genes tested. The distribution of gNAS scores (A,B) and gFib scores (C,D). Plots (A,C) show the distribution of gNAS or gFib scores for the top 1000 genes in each sample. Plots (B,D) show the relationship between mean gNAS and gFib scores and histological assessments.

Figure 4

Lasso regression of gene expression values against mean gNAS (A,C,E) or gFib scores (B,D,F). Figures (A,B) show the results of 5-fold cross-validation for each model, which have 19 and 18 predictors, respectively. The strong performance of the models in cross-validation demonstrates that disease severity can be assessed from the expression levels of a relatively small number of genes. Figures (C,D) provide the scaled variable importance for model predictors. Figures (E,F) show the standardized regression coefficients for each model.

Figure 5

Patterns of gNAS and gFib scores across patient samples reveal distinct molecular profiles. Panels (A,B), respectively, show standardized gNAS and gFib scores across sets of genes that were identified through gene shaving. Sample clusters in these panels show distinct patterns regulation across these genes, and thus represent patients with distinct molecular profiles. Panels (C,D) show the distributions of mean standardized scores for each sample cluster. Within these plots, patterns across gene clusters (x-axis) represent the average molecular profiles of the sample clusters. Panel (E) shows the intersection of the gNAS- and gFib-based sample clusters and provides the number of samples in each cluster pair. Simultaneous consideration of the two partitions provides additional granularity in sample classification. Panels (F,G) show the most strongly represented Reactome pathways in each gene cluster (by Fisher’s exact test). The pathways represented are closely linked to NAFLD progression.

Figure 6

Patterns of pathway-level regulation with respect to gNAS and gFib scores. The heatmaps (A,B) shows the clustering pattern of samples (columns) and MSigDB hallmark pathways (rows) with respect to mean gene-level scores (values represent column-wise Z-scores). Sample clusters show distinct patterns of pathway-level regulation. Panels (C,D) show the mean sample-wise Z-score in each cluster for the gNAS and gFib analyses, respectively. Higher values in both figures are consistent with relatively advanced disease states. Panel E shows the intersection of the gNAS- and gFib-based sample clusters and provides the number of samples in each cluster pair.