Nonalcoholic fatty liver disease (NAFLD) is the commonest cause of chronic liver disease worldwide and a growing healthcare burden. The pathobiology of NAFLD is complex, disease progression is variable and unpredictable, and there are no qualified prognostic biomarkers or licensed pharmacotherapies that can improve clinical outcomes; it represents an unmet precision medicine challenge. We established a retrospective multicentre national cohort of 940 patients, across the complete NAFLD spectrum, integrating quantitative digital pathology, hepatic RNA-sequencing and 5.67 million days of longitudinal electronic health record follow-up into a secure, searchable, open resource (SteatoSITE) to inform rational biomarker and drug development and facilitate personalised medicine approaches for NAFLD. A complementary web-based gene browser was also developed. Here, our initial analysis uncovers disease stage-specific gene expression signatures, pathogenic hepatic cell subpopulations and master regulator networks associated with disease progression in NAFLD. Additionally, we construct novel transcriptional risk prediction tools for the development of future hepatic decompensation events.
Resource
SteatoSITE: an Integrated Gene-to-Outcome Data Commons for Precision Medicine Research in NAFLD
https://doi.org/10.21203/rs.3.rs-2805134/v1
This work is licensed under a CC BY 4.0 License
Journal Publication
published 30 Oct, 2023
Version 1
posted
You are reading this latest preprint version
Nonalcoholic fatty liver disease (NAFLD) is a growing global public health concern, affecting more than one in four adults (over two billion people) and imposing a substantial burden of ill health and socioeconomic problems1,2. It is more common in patients with obesity and type 2 diabetes, and is associated with increased morbidity and mortality from cardiovascular disease, chronic kidney disease and many cancers. NAFLD is an umbrella term used for a range of conditions caused by a build-up of fat in the liver. This includes isolated fatty liver (simple steatosis) and nonalcoholic steatohepatitis (NASH), the inflammatory subphenotype of NAFLD that can lead to progressive liver fibrosis, cirrhosis, and hepatocellular carcinoma (HCC). Notably, over the past ten years, there has been a tenfold increase in the number of patients with NAFLD-related cirrhosis who require a liver transplant3 and recent Bayesian modelling forecasts a staggering global prevalence rate of 55.7% by 20404.
The pathogenesis of NAFLD is complex and multifactorial5,6. The primary insult of excess lipid accumulation (“substrate overload”) is followed by variable contributions from pathogenic drivers including lipotoxicity and immune-mediated inflammation, and disease progression further enhanced/suppressed by modifiers such as genetic variants7 and gut microbiota dysbiosis8,9. Crucially, not everyone with a fatty liver goes on to develop adverse liver-related events or death, but there is a critical lack of prognostic biomarkers to enable innovations in clinical care pathways and personalised approaches in NAFLD. The importance of histological fibrosis in predicting outcomes in NAFLD has been highlighted by several studies10–12, with discrete fibrosis stages dramatically influencing the risk of all-cause mortality and liver-related morbidity and mortality. Yet fibrosis progression rates and outcome prediction in patients remain uncertain and clinical trials that have been used post hoc to understand the natural history of NAFLD are limited by selection bias and rigid entry criteria13.
Numerous drugs with plausible biological targets have failed to show robust efficacy in late-phase interventional trials in NASH, especially in patients with cirrhosis14, adding to the so-called “NASH graveyard”. Consequently, the field is shifting focus to combination drug regimens and, increasingly, researchers are leveraging patient-centred big data to bolster drug discovery efforts15. However, the most efficacious drugs/combinations to mitigate adverse clinical outcomes in NAFLD remain to be defined.
To address these issues and provide a resource suitable for the development of precision medicine strategies in NAFLD, we established SteatoSITE – the world’s first data commons16 for NAFLD research. SteatoSITE integrates multimodal multiscale retrospective data, at a national level, from 940 cases across Scotland consisting of: quantitative histopathological assessment of archival liver tissues, hepatic bulk RNA-sequencing, and routine clinical data extracted from electronic health records (EHR) and other administrative sources. Distinctively, SteatoSITE is outcome-rich, with data curated in a longitudinal time-stamped format comprising 5.67 million days of clinical follow-up, enabling translational analyses to elucidate the pathogenesis of NAFLD and accelerating the discovery of potential prognostic biomarkers and tractable therapeutic targets. Here, we describe and utilise this unique resource by providing a comprehensive analysis of hepatic gene expression and cell type composition across the complete NAFLD spectrum, by identifying transcription factor (TF)-regulated gene networks (regulons) associated with disease progression, and by developing a novel transcriptome-based risk prediction model for hepatic decompensation.
Cohort composition
A total of 940 cases from the three participating NHS Scotland Biorepositories (Lothian, Greater Glasgow & Clyde, and Grampian) were included, representing the full histological spectrum from normal liver tissue to NAFLD-related cirrhosis (Fig. 1a). Demographic and phenotypic characteristics of the cohort are shown in Table 1. Cases with a liver tissue sample acquired between January 2000 and October 2019 were selected. All patients were ≥ 18 years of age at the tissue sampling date. Data from EHRs and national datasets were retrieved, where available, from a period between ten years before the tissue sampling date until May 2020, comprising over 5.67 million days (or approximately 15,547 years) of comprehensive routine clinical data and annotated timelines created (Fig. 1b, Extended Data Fig. 1, Supplementary Table 1). The post-tissue sampling periods for biopsies, resections, and explants are shown in Extended Data Fig. 2a.
Histopathological characterisation of the SteatoSITE cohort
Traditional subjective histopathology assessments
Scans of haematoxylin and eosin (H&E) and picrosirius red (PSR)-stained sections from each case were assessed by one of three consultant pathologists with a specialist interest in liver pathology. From the H&E-stained sections, NAFLD activity scores (NAS), components of the Steatosis, Activity and Fibrosis (SAF) scoring system, and a score-independent clinical histological diagnosis of NASH were given. From the PSR-stained sections, a NASH-Clinical Research Network (NASH-CRN) and modified Ishak score for staging of fibrosis were given. The inter-rater agreement levels of the scoring pathologists (Supplementary Table 2), assessed on a set of 20 cases after the scoring harmonisation exercise, were at the upper end of the expected ranges17,18. Spearman’s rho correlation coefficient was used to assess the relationship between modified Ishak fibrosis scores and NASH-CRN fibrosis scores, demonstrating the expected significant positive correlation (rs = 0.98, P < 2.2e-16).
Of the 940 cases, 659 were needle biopsies, 227 were from hepatic resections, and 54 were explants for end-stage NAFLD cirrhosis. A clinical histological diagnosis of NASH was given in 455 (48.4%) of cases (Extended Data Fig. 2b). Approximately equal numbers of cases were scored at each point of the NASH-CRN fibrosis scale. The relationship of the NAS components and the NASH-CRN fibrosis score are shown in Fig. 1c.
Quantitative digital pathology feature measurement
A pixel classifier was trained to identify fat and PSR-positive tissue within scans of PSR-stained sections (Extended Data Fig. 2c). All classified images were manually quality-controlled and any images containing large fragments of liver capsule, large portal tracts or hilar tissue, or with artefact that was easily ignored during subjective assessment but erroneously computationally classified were removed. The derived fat percentage of the tissue correlated with the assigned steatosis score of the NAS (and SAF) systems, as expected (Spearman’s rho correlation coefficient rs = 0.53, P < 2.2e-16, Extended Data Fig. 2d). The derived PSR-positive percentage of the tissue also correlated with both NASH-CRN (Spearman’s rho correlation coefficient rs = 0.61, P < 2.2e-16, Extended Data Fig. 2e) and modified Ishak (Spearman’s rho correlation coefficient rs = 0.61, P < 2.2e-16, Extended Data Fig. 2f) assigned scores.
Association of fibrosis stage with all-cause mortality and liver-related outcomes
The extensive annotation of individual case timelines with clinical data, anchored by the specimen date and retrieved from the EHR and national datasets, allowed the relationship between data derived from the histological sections and patient outcomes to be examined. Using the 659 needle biopsy cases only, the clear and expected relationship between assigned NASH-CRN fibrosis stage and all-cause mortality was observed (Fig. 1d). Stepwise increases in mortality were evident through progression from stages F0 to F4. An unbiased algorithmic approach to cluster survival curves19 created two clusters (F0,1,2) and (F3,4), supporting the hypothesis that fibrous bridging of vascular structures is a critical pathophysiological event with prognostic importance in progressive fibrogenesis.
Standard Cox regression modelling of all-cause mortality using age at biopsy, sex, and NASH-CRN fibrosis stage as covariates showed that age and NASH-CRN fibrosis stage F4 were independently associated with higher all-cause mortality (Supplementary Table 3).
SteatoSITE is a secondary care, outcome-enriched cohort. The annotated patient timeline of 291 cases, including 183 of 659 biopsy cases, contained an outcome coding for at least one of the events defined by expert consensus to represent decompensation in cirrhotic patients20 or using UK Operations/Procedure (OPCS4) coding data identifying activity relating to cirrhosis-related hospital admissions21. This high number of events contrasts with a recently published prospective observational study12. Using the 106 cases where an event coding related to decompensation was not present on the patient timeline prior to the biopsy, F3 and F4 NASH-CRN fibrosis stage was predictive of subsequent decompensation (Fig. 1e). Kaplan-Meier estimator curves for liver-related events (excluding death) associated with F0, F1, and F2 were placed in the same cluster by an unbiased clustering approach, but those associated with F3 and F4 were distinct.
Cox regression modelling of hepatic decompensation events on the cause-specific hazard, with death as a competing risk, using age at biopsy, sex, and NASH-CRN fibrosis stage showed that NASH-CRN fibrosis stages F3 and F4 were associated with increased decompensation events (Supplementary Table 4); Fine-Gray regression for the proportional hazards modelling of the subdistribution hazard is also reported. Median intervals to decompensation or censoring (death or end of follow-up) of included cases are shown in Supplementary Table 5.
Finally, the development of HCC is a low frequency outcome in NAFLD. The complete SteatoSITE cohort includes 80 patients with a coding of HCC at any point in their annotated timeline. The increased risk of HCC in patients with non-cirrhotic NAFLD is also well-recognised22,23. The SteatoSITE cohort includes 227 resection specimens, with HCC being present in the annotated patient timeline in 44 of these. Notably, 36 of these cases did not have histopathological evidence of cirrhosis in the background liver (Extended Data Fig. 2g). Focusing on the liver biopsy cases where no event coding for HCC was present prior to the biopsy, ten patients received a subsequent coding of HCC within the follow-up period included in the data commons, and assigned NASH-CRN F4 fibrosis stage at biopsy was predictive of the subsequent development of HCC (Fig. 1f).
Overall, these data clearly demonstrate that within the SteatoSITE cohort, liver fibrosis stage is predictive of future all-cause mortality, hepatic decompensation events and HCC development.
Transcriptomic profiling across the NAFLD spectrum uncovers gene signatures for steatohepatitis and fibrosis
The SteatoSITE cohort includes high-quality hepatic RNA-seq data in 707 out of the 940 total cases (comprising 538 biopsies, 39 explants and 130 resections), after applying pre-specified quality control criteria appropriate for archival FFPE samples, including DV200 > 30%24. Overall, larger liver resection tissues yielded poorer RNA quality compared with much smaller needle biopsy specimens, likely related to sample fixation.
Normal liver controls (n = 34) were also included in SteatoSITE for comparative RNA-seq analysis. To confirm the suitability of these control samples, we showed that their transcriptional profile strongly correlated with normal liver samples from independent hepatic RNA-seq datasets25,26 (average Spearman’s correlation rs > 0.9, adjusted P < 2 e-16; Extended Data Fig. 3).
We performed RNA-seq variant calling to detect single nucleotide polymorphisms (SNPs) previously identified as risk modifiers of NAFLD progression (“rs738409”, “rs72613567”, “rs58542926” and “rs641738” in genes PNPLA3, HSD17B13, TM6SF2 and MBOAT7, respectively). Sequencing coverage of these genes, and therefore the ability to call variants, was variable (Supplementary Table 6). However, the rs738409 C > G p.I148M variant in PNPLA3 (the strongest genetic risk factor for NAFLD and its severity) was successfully called in 646 of 707 samples with RNA sequencing available. The prevalence of GG, GC, and CC genotypes amongst cases where a call was possible was 17.02%, 26.63% and 55.42%, respectively. After controlling for fibrosis stage, genotype status for pre-defined variants had no significant effect on outcomes (data not shown).
For differential gene expression analysis, samples were grouped according to NAS and different fibrosis stages (F0/F1, F2, F3 and F4) and compared with normal liver control samples (Extended Data Fig. 4). The data structure was examined by principal component analysis (PCA). As shown in Fig. 2a, the first two PCs explained 8.86% and 7.94% of the observed variation in gene expression. There was modest segregation according to the fibrosis stage.
When samples were classified according to fibrosis stage, we found 218, 478, 1114 and 2800 differentially expressed protein-coding genes when comparing F0/F1, F2, F3 and F4 liver with normal liver controls, respectively, and 99 genes differentially expressed in common across stages versus controls (Fig. 2b). Gene set enrichment analysis (GSEA) of the upregulated hepatic genes showed that enriched gene ontology (GO) terms associated with early fibrosis stages included cytokine activity and structural components (Extended Data Fig. 5). A representation in two-dimensional semantic space of the GO terms mapped to those genes differentially expressed in F4 versus control, with selected terms highlighted, is shown in Fig. 2c; this includes GO terms related to the ‘immune system’ in addition to those related to ‘wound healing’. The most differentially expressed genes mapped to GO terms ‘collagen fibril organization’, ‘extracellular matrix assembly’, ‘wound healing’, and ‘vascular wound healing’ are shown in Fig. 2d.
Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways linked to extracellular matrix were enriched from stage F2 onwards, as expected. Notably, other upregulated pathways were those linked to cardiomyopathies (q < 0.05; Fig. 2e). In contrast, downregulated genes were enriched for pathways associated with signalling, predominantly in early disease (F0/1), and fatty acid metabolism and peroxisome pathways were enriched from F2 to F4 (q < 0.05; Fig. 2f).
When the samples were grouped according to NAS, we identified 203 DEGs when comparing NAS ≤ 2 (“low NAS”) with controls, 621 in the NAS 3–4 group (“borderline NAS”) and 793 in the NAS > 4 group (“high NAS”); 182 genes were shared across all NAS groups (Extended Data Fig. 6a). GSEA revealed, among others, significantly enriched GO terms that participate in ‘translation’ and the ‘immune system’ (q < 0.05; Extended Data Fig. 6b-d).
Hence, the SteatoSITE data commons enables a comparison of whole liver gene expression profiles according to histological disease stage. To maximise the accessibility and utility of this resource, we developed an open-access gene browser (https://shiny.igc.ed.ac.uk/SteatoSITE_gene_explorer/), allowing high-level assessment and visualisations of user-selected gene expression according to fibrosis stage.
Characterisation and predictive value of cell type compositions from SteatoSITE bulk RNA-sequencing data using single-cell deconvolution
Using single-cell RNA-seq (scRNA-seq), we have previously identified distinct populations of scar-associated monocyte-derived macrophages, mesenchymal cells and endothelial cells which populate the fibrotic niche in patients with advanced cirrhosis and interact to regulate liver fibrogenesis27. However, the presence of specific pathogenic cellular subpopulations across the full NAFLD disease spectrum and how these cells relate to clinical outcomes remains uncertain. To address this, we performed deconvolution of the bulk RNA-seq SteatoSITE data using a publicly available fully annotated scRNA-seq reference dataset compiled from healthy and cirrhotic patients27 to estimate the proportions of specific hepatic cell types in each SteatoSITE sample and to correlate these with histopathological features and patient outcomes.
This analysis demonstrated that the proportion of hepatic scar-associated macrophages (SAMac), a key regulator of liver fibrosis27, correlated positively with liver fibrosis stage and NASH activity across the full NAFLD spectrum (Fig. 3a). In contrast, tissue resident macrophage (Kupffer cell) proportions declined significantly with more advanced fibrosis (Fig. 3a). We validated these transcriptomic findings at protein level, using a bespoke MultiOmyx liver multiplex immunofluorescence (mIF) assay in an independent NAFLD cohort (n = 43), confirming a significant positive correlation of hepatic HLA-DR+CD9+CD14+ SAMac numbers with liver fibrosis stage (Fig. 3b). Positive correlations between mIF SAMac numbers and NASH steatosis, lobular inflammation and ballooning scores were less pronounced than fibrosis stage (Extended Data Fig. 7a-c), mirroring our findings from the deconvolution analysis (Fig. 3a).
Deconvolution analysis also revealed significant positive correlations between fibrosis and the proportion of scar-associated mesenchymal cells, hepatic arterial endothelial cells, and lymphatic endothelial cells (Fig. 3a), in keeping with the key roles for mesenchymal cell activation and progressive arterialisation of the hepatic microcirculation with loss of normal specialised liver sinusoidal endothelial cell (LSEC) phenotype (“capillarisation”) in the development of liver fibrosis across the full NAFLD spectrum. Interestingly, expansion of plasma cells was also associated with more advanced hepatic fibrosis (Fig. 3a).
Finally, the proportions of cell types, derived by single-cell deconvolution analysis of SteatoSITE bulk RNA-seq data, were examined for their value in predicting adverse clinical outcomes. Strikingly, increased hepatic proportions of SAMac, scar-associated mesenchymal cells, hepatic arterial endothelial cells, lymphatic endothelial and plasma cells were predictive of higher future all-cause mortality (Fig. 3c) and hepatic decompensation events (Fig. 3d). In contrast, higher proportions of more homeostatic liver resident cell types such as vascular smooth muscle cells and LSECs was protective against future mortality or hepatic decompensation (Fig. 3c-d). Hence, in addition to histology and bulk transcriptomics, changes in the cellular composition of the liver may offer key prognostic information in patients with NAFLD.
Construction of a transcriptional risk prediction model for hepatic decompensation events
The annotated patient timelines in SteatoSITE allow predictive tools to be developed. To demonstrate, we used the SteatoSITE RNA-seq data and associated clinical outcomes to develop a novel transcriptome-based risk prediction model for hepatic decompensation. Such transcriptional risk scores (TRSs) based on transcript abundance are physiologically closer to the phenotype of interest, require smaller training samples and offer greater portability across diverse ancestry groups than polygenic risk scores (PRSs) using genomic variants (single-nucleotide polymorphisms (SNPs))28,29.
As the initial feature set, we used the 1127 protein-encoding genes that were differentially expressed in advanced (F3/F4) compared with early (F0/F1) stage disease (P < 0.05, log fold change (FC) ≥ 1). Univariate Cox regression identified 955 DEGs as significantly related to decompensation events (P < 0.01). To develop a sparse feature set, 10 runs of a 10-fold LASSO regularised Cox regression were performed and the coefficients for selected genes were derived (Extended Data Fig. 8a-b). The final TRS predicting hepatic decompensation was composed of the expression of 15 genes: metallothionein-1F (MT1F), potassium voltage-gated channel subfamily H member 7 (KCNH7), collagen alpha-1(XXV) chain (COL25A1), GTP-binding protein Rhes (RASD2), connective tissue growth factor (CTGF), stanniocalcin-1 (STC1), glial cell line-derived neurotrophic factor (GDNF), paired related homeobox 1 (PRRX1), fibroblast growth factor 7 (FGF7), lipocalin like 1 (LCNL1), dipeptidase 1 (DPEP1), chordin like 2 (CHRDL2), LIM homeobox protein 6 (LHX6), POU domain class 4 transcription factor 1 (POU4F1), cadherin 16 (CDH16) (Expression across fibrosis stages as plotted by the SteatoSITE gene browser shown in Extended Data Fig. 9).
The risk scores were calculated using the formula indicated in the Methods. According to the median risk score, patients were split into high- and low-risk groups (Fig. 4a-b, Extended Data Fig. 8c). Interestingly, samples with a higher risk score not only had a higher fibrosis stage but also a higher hepatocyte ballooning score (Fig. 4c). Time-dependent ROC curves were used to assess the predictive ability of the model at specific times after biopsy; the AUROCs were 86.24% (SE 5.11), 80.97% (SE 4.62) and 83.26% (SE 3.86) for 1-, 3- and 5-year risk of decompensation events, respectively (Fig. 4d).
All the biopsy samples could be stratified using the TRS into high- and low-risk groups with significantly different decompensation trajectories. Over post-biopsy follow-up of up to 20 years, those with a high-risk TRS had a cumulative decompensation event probability of over 0.65 compared with a cumulative probability of 0.11 in those with a low-risk TRS (Fig. 4e). Next, to derive additional prognostic information beyond routine fibrosis staging, samples with mild (F0/F1) or advanced (F3/F4) scarring were stratified using the TRS into groups at high- or low-risk of future decompensation (Fig. 4f). Although there were insufficient F3/F4 patients with a low TRS to enable further categorisation of these patients, application of the TRS augmented risk stratification in NASH patients with early-stage fibrosis.
Master regulator analysis identifies thyroid receptor hormone beta as a critical suppressor of disease progression
We sought to further exploit the rich RNA-seq dataset using the high-risk genes identified as prognostically important components of the TRS to derive additional biological understanding of the related transcriptional regulatory network (TRN) in NAFLD. The TRN, consisting of TFs and regulated target genes, was inferred from the whole gene expression dataset. Regulons, a set of genes regulated by a specific TF, were constructed for all TFs catalogued by Lambert et al30. A regulon activity score for each sample was estimated using a two-tailed GSEA approach. Regulons cluster based on activity broadly into two groups, those with high activity in advanced fibrosis and low activity in early stages, and those showing the opposite pattern (Fig. 5a).
To identify which transcriptional networks lay upstream of the expression of high-risk genes representing the TRS, master regulator analysis31 was undertaken to identify the statistical significance of the overlap between the regulon (gene targets) of each TF and the 15 genes of the TRS, corrected for multiple hypothesis testing. Three regulons (gene networks regulated by adipocyte enhancer-binding protein 1 (AEBP1), thyroid hormone receptor beta (THRB), and zinc finger protein basonuclin-2 (BNC2)) contained significantly greater numbers of the TRS genes than expected by chance, suggesting that these three networks may be critical to NAFLD disease progression leading to decompensation events.
To examine the patterns of activity of these regulons during progression of disease stage, the mean activity scores for each regulon from each NASH-CRN fibrosis stage were used as pseudo-timeseries data for unsupervised soft clustering. AEBP2 regulon activity is in a cluster of regulons whose activity is low in F0-2 but increased in F3 and F4 whereas BNC2 regulon activity increases more uniformly from F0 to F4. In contrast, THRB regulon activity is high in F0 and F1 but decreases as scarring progresses from F2 to F4 (Fig. 5b).
The relationship between the inferred AEBP1, THRB, and BNC2 regulons is shown in Fig. 5c, annotated with the direction of regulation of each gene target. Most gene targets for AEBP1 and BNC2 are positively regulated whilst THRB exerts largely negative regulation of gene targets shared with both AEBP1 and BNC2.
Given the potential suppressive role for THRB over disease-promoting core gene targets, and because THRB agonists are in current clinical trials for NASH, we examined THRB regulon activity in more detail. In biopsy cases with no hepatic decompensation prior to biopsy, differential regulon activity was estimated, and the predictive value of THRB regulon activity for decompensation events examined (Fig. 5d). A ranked THRB regulon differential activity plot confirms the negative relationship between THRB regulon activity and disease stage in the biopsy-restricted subset and the Kaplan-Meier plot indicates that lower THRB regulon activity is predictive of hepatic decompensation.
However, the predictive value of THRB regulon activity may merely be through association with disease stage. To determine if THRB regulon activity was predictive of hepatic decompensation within samples of matched fibrosis stages, maximally selected rank statistics were used to determine the optimal THRB regulon activity cutpoint that divided samples into high- and low-risk groups. Using this cutpoint on histologically low-risk samples with F0 and F1 fibrosis, THRB regulon activity identified a subset of NAFLD patients with high risk of hepatic decompensation (Extended Data Fig. 10). Deriving a new cutpoint within only the histologically advanced stages F3 and F4, allowed stratification into two distinct groups with either rapid or slower progression to a decompensation event (Fig. 5e).
The potential for a personalised approach, for clinical trials or routine practice, using estimated THRB regulon activity is illustrated by case studies from the SteatoSITE dataset. Two cases scored F4 (cirrhosis) on PSR-stained sections with high and lower estimated THRB activity are shown in Fig. 5f; the patient with estimated THRB regulon activity below the derived threshold (top) had a coded decompensation event 224 days after biopsy but the patient with THRB activity above the cutpoint (bottom) did not experience a coded decompensation event in over 4500 days of follow-up.
At the other end of the disease stage spectrum, two cases that scored F0 on the PSR-stained sections are shown in Extended Data Fig. 10; the case with high estimated THRB regulon activity derived from RNA-seq (bottom) was not associated with hepatic decompensation whereas the case with low estimated THRB regulon activity (top) was.
We have established the world’s first data commons for NAFLD research. This unique human resource - comprising quantitative histology, transcriptomics, and rich clinical data across the full disease severity spectrum - is a dataset of broad utility that will enable discovery science, digital pathology research, artificial intelligence approaches and translational studies in NAFLD. SteatoSITE is complementary to ongoing global initiatives in NAFLD such as biomarker-focused research being conducted by the LITMUS and NIMBLE consortia. However, crucially, gene signatures and histopathological assessments in SteatoSITE are linked to clinical outcomes to enable rational diagnostic and treatment approaches for patients with different stages and activity of disease to be determined. Indeed, a recent authoritative expert review stated that “for more refined treatment of NASH, orthogonal approaches that integrate genetic, clinical, and pathological data sets may yield treatments for specific subphenotypes of the disease”24. SteatoSITE provides the relevant multiscale and multimodal data required to meet this challenge.
Consistent with previously published observational data32, we began by showing the strong positive association between histological fibrosis stage and future clinical outcomes including all-cause mortality, HCC, and a composite outcome of hepatic decompensation events. These findings validated the broad content and data integration of the SteatoSITE cohort. We went on to perform a comprehensive bioinformatic analysis of the hepatic bulk RNA-seq data, identifying key genes and enriched gene functions/pathways that characterise discrete histological subtypes and generating disease stage-specific transcriptional profiles to inform drug target discovery and the development of clinically relevant biomarkers in NAFLD. In addition, we provide an open-access web-based gene browser so the scientific community can interactively explore and visualise the SteatoSITE RNA-seq data, thus increasing the accessibility and impact of our work.
The advent of single-cell transcriptomic technologies is transforming our understanding of the pathobiology of liver diseases including NAFLD, identifying key pathogenic cell types and candidate therapeutic targets33–35. However, the study of these cell populations across the full disease spectrum remains limited, with minimal data defining which specific cell states are associated with adverse clinical outcomes. Here, we employed a single-cell deconvolution approach to interrogate changes in cellular composition in the SteatoSITE cohort, confirming expansion of SAMac, scar-associated mesenchymal cells, plasma cells, hepatic arterial endothelial cells, and lymphatic endothelial cells during NAFLD progression. Crucially, expansion of these populations also correlates with poor long-term patient prognosis, highlighting their potential role as therapeutic targets in NAFLD and informing novel immunohistochemical approaches to better characterise liver biopsy tissue and identify high-risk NASH subphenotypes. Indeed, multiplex immunohistochemistry is already being utilised to improve stratification in patients with a range of cancers36–38, so future refinement of our liver macrophage-focussed high-plex MultiOmyx assay to include specific markers for pathogenic mesenchymal, endothelial and plasma cells will enable a comprehensive assessment of the liver fibrotic niche from a single slide and may facilitate improved identification of high-risk NAFLD patients than existing histological approaches.
As the individual patient-level datapoints in SteatoSITE are temporally defined, researchers can perform time-to-event predictions using survival analysis methods. Here, we illustrated this by showing that hepatic gene expression data in our cohort had a high predictive value for risk of future hepatic decompensation events, over and above standard fibrosis staging. The transition from compensated cirrhosis to decompensated cirrhosis occurs at a rate of about 5–7% per year39. Decompensation represents a key prognostic inflection point in the natural history of chronic liver disease, as the median survival drops from more than 12 years for compensated cirrhosis to about 2 years for decompensated cirrhosis39. However, in patients with advanced fibrosis or cirrhosis (F3/F4), individual risk of decompensation is variable and hard to predict, and preventative therapies are lacking. Based on automated variable selection methods to reduce overfitting, we identified a 15-gene panel and derived a risk score which could classify stage-defined NAFLD patients into high- or low-risk groups for decompensation40,41. To the best of our knowledge there are no TRS reported in NAFLD but, notably, a TRS approach outperformed static genotypic risk assessment in distinguishing Crohn’s disease from healthy samples and predicting complications42. Genes represented in our TRS include five members predicted to be secreted proteins according to The Human Protein Atlas (https://www.proteinatlas.org), which could be explored further in suitably well-annotated clinical samples. Although not previously reported in NAFLD, CHRDL2 is upregulated in a range of tumour tissues including HCC43, while STC1, CTGF44, GDNF and FGF7 are all linked to hepatic fibrogenesis45,46. Additionally, STC1 is a potential serum biomarker for hepatitis B virus-associated liver fibrosis45 and secreted FGF747 may be a prognostic marker in cholangiocarcinoma. The predictive model requires external validation, ideally in prospective longitudinal studies. However, to our knowledge, there are no accessible NAFLD datasets with integrated pathology, transcriptomics, and abundant clinical outcome data to enable this.
Finally, to search for new and highly effective upstream therapeutic targets, we used high-risk genes from the TRS to uncover core gene networks and master regulators likely to exert influence over disease progression and patient outcomes. Intriguingly, one regulon – THRB – was suppressive of the other two identified. Moreover, THRB regulon activity not only decreased with advancing fibrosis stage but also predicted future hepatic decompensation (beyond standard fibrosis scoring). Liver-targeted thyromimetics selectively activating THRB such as Resmetirom (MGL-3196) and VK2809 increase hepatic fat metabolism and reduce lipotoxicity and have recently emerged as leading pharmacological candidates for the treatment of NASH48. Indeed, Resmetirom49 has advanced to phase 3 trials and a NASH cirrhosis outcomes trial is underway. Our data reinforce the importance of THRB as a therapeutic target with the potential to impact hard clinical endpoints and illustrate how adding transcriptomic information to histology might further stratify prognosis and enable tailoring of treatment.
The creation of SteatoSITE was facilitated by the Scottish demographics and healthcare infrastructure. Scotland has a well-established ecosystem for precision medicine consisting of a stable population base of ~ 5.5 million people, a significant incidence of major chronic disease (death rates from chronic liver disease in Scotland are 70% higher than the UK average50), a single healthcare provider (National Health Service (NHS) Scotland) and an EHR system operating on a national scale with a unique common field identifier (Community Health Index (CHI) number) allowing characterisation and longitudinal follow up of well-defined patient cohorts; these key elements differentiate Scotland from many other countries.
Given the scale and complexity of SteatoSITE and its use of routine multicentric clinical data, several challenges and limitations are acknowledged. Leveraging EHR data for NAFLD research is recognised as potentially powerful, but consensus guidelines have only recently emerged20,51. Using three regional NHS Safe Havens, we encountered several problems that required an exacting and systematic data cleaning process, including: data duplication, human error (e.g., incorrect data entry, typographical errors, sample mislabelling), issues with data standardisation, errors due to different delimiters/encoding in input files, data formatting, and missing data/completeness. Notably, some source data (e.g. Fibroscan results) are inconsistently recorded across health boards and were therefore hard to capture, even using laborious manual methods. Additionally, there are specific limitations regarding the cohort itself. Inevitably, using a secondary care tissue-first selection process introduces inherent spectrum bias. This is a strength in terms of outcome enrichment but means that SteatoSITE will have less value for modelling the population-level natural history of NAFLD. The demographic make-up of the cohort also lacks ethnic diversity (majority white Scottish), so caution is advised regarding the generalisability of findings to other geographical areas and ethnic populations. Finally, confounding due to unsuspected52 or poorly recorded alcohol use is a perennial issue in NAFLD research. We were fastidious in our case selection and collected ICD codes that can identify alcohol use disorders recorded in the EHR20. Outside of a rigidly controlled clinical trial, the real-world contribution of mild-to-moderate alcohol consumption in observational studies is difficult to eliminate, whilst its impact on NAFLD also remains controversial53,54.
SteatoSITE has laid the groundwork for a variety of important research themes in NAFLD, providing a rich multimodal database to inform drug and biomarker development and a catalogue of high-resolution whole-slide images as a valuable testbed for novel digital pathology algorithms. Here we describe SteatoSITE version 1, but our intention is that this data commons will be continuously added to and refined to enhance its content and operability.
Regulatory approvals
Anonymised tissue was supplied after approval by the National Health Service Research Scotland (NRS) Biorepository network. Unified transparent approval for data inclusion in this pan-Scotland project was provided by the West of Scotland Research Ethics Committee 4 (Reference: 20/WS/0002; 18th February 2020), Public Benefit and Privacy Panel for Health and Social Care (PBPP; Reference: 1819-0091; 4th June 2021), Institutional Research & Development departments and Caldicott Guardians.
Study cohort and sample collection
Initial case selection was based on the availability of archival liver tissue (from biopsies, resections, or explants that were surplus to diagnosis) in formalin-fixed paraffin-embedded (FFPE) blocks available within the NRS Biorepository network, with the clinical and/or histological diagnosis of NAFLD, and meeting the following inclusion/exclusion criteria:
Inclusion criteria - male or female, > 18 years of age at the time of tissue sampling, all ethnic groups, socio-economic backgrounds and health status, dead or alive at the time of inclusion into data commons.
Exclusion criteria - cases were excluded if any of the following applied: documented history of chronic liver disease of any non-NAFLD aetiology, including alcohol-related liver disease, chronic viral hepatitis, haemochromatosis, Wilson’s disease, autoimmune hepatitis, primary biliary cholangitis, primary sclerosing cholangitis, and patients with excessive alcohol use documented within the clinical data supplied on the specimen request form (> 21 units/week for men, > 14 units/week for women) or histological features suggesting a secondary non-NAFLD diagnosis.
FFPE blocks of normal non-lesional liver from 40 control cases, defined as liver samples without fat deposition or any active primary parenchymal disease, were identified from resection specimens.
Clinical data extraction and processing
After approval from the NHS Scotland PBPP, we obtained relevant unconsented clinical data that included: demographics (age, gender, ethnicity, Scottish Index of Multiple Deprivation (SIMD)), diagnostic coding (International Classification of Diseases (ICD)-9, ICD-10) for inpatient and outpatient episodes and procedures (OPCS Classification of Interventions and Procedures (OPCS-4)), routine laboratory blood tests (biochemistry and haematology), liver transient elastography (Fibroscan®) where available, prescribing information, cancer and death registry data. The comprehensive list of variables (including ICD codes) is shown in Supplementary Table 1. To maximise future comparability and generalisability of results across studies, we followed recent expert consensus guidelines for using administrative coding in EHR-based research of NAFLD20.
The EHR data was uploaded to the Precision Medicine Scotland-Innovation Centre Healthcare Landing Zone from three NHS Safe Havens as comma separated value (.csv) files. In areas where the Safe Havens had poor coverage from their national or regional datasets, additional data files were provided by the Biorepository network. In total, 38 csv files containing 1,450 columns and 1,055,618 rows were supplied for processing. Upon upload, data was quality controlled, processed and further minimised to comply with regulatory requirements and ethical approvals.
A Python script was developed to carry out these data processing tasks. These included data cleaning, data restructuring and the standardisation of file inputs to resolve differences in file encoding, text delimitation and formatting. Where appropriate, differences in data labelling were resolved by adopting the naming conventions used by the Glasgow Safe Haven. After standardising case and spacing, 962 distinct columns remained from the original 1,450 columns. After a more detailed review, columns containing overlapping clinical information were relabelled and merged. This further reduced the dataset to 534 distinct columns which are present in the final project database.
A final key requirement of the pseudonymisation processing was to ensure that all event dates were obfuscated to minimise re-identification risk whilst still allowing for clinically meaningful patient timelines to be established. To achieve this, the duration of time between the specimen acquisition date and the clinical event recorded in the EHR was calculated as Δt in seconds and stored.
Histopathology methods
Tissue preparation
Blocks for which sufficient tissue was available were assigned a Study ID, and samples were cut for histology (3 x 4 µm), and where possible, RNA extraction (4 x 10 µm). Unstained sections were sent for histological staining by the NHS Lothian pathology laboratory and curls for RNA extraction at the Genetics Core in the Edinburgh Clinical Research Facility. Only NRS biorepository staff independent from the study group held the participant key; they passed data linkage documents to the local NHS Safe Havens who were also independent of the study group. No information from the original clinical pathology report was supplied with the sections.
One section from each case was haematoxylin and eosin (H&E)-stained and one section was picrosirius red (PSR)-stained. Stained sections of control cases were reviewed to ensure no histopathological abnormalities were present. Stained sections of NAFLD-spectrum cases were scanned at x20 depth on a Hamamatsu NanoZoomer whole-slide scanner.
Histopathological scoring
Raw .ndpi files of H&E and PSR-stained NAFLD-spectrum cases were uploaded to a custom server installation of OMERO55 hosted by PMS-IC (Glasgow, UK).
In June 2019, prior to SteatoSITE case scoring, three consultant pathologists and participants in the UK National Liver External Quality Assurance scheme (TJK, GK, PK) undertook an in-person scoring harmonisation exercise where examples of histological features scored in the NAS17 and SAF56 scoring systems and for staging using NASH-CRN and modified Ishak57 systems were reviewed. To assess inter-rater performance following this meeting, scans of H&E- and PSR-stained sections from 20 NAFLD-spectrum cases including biopsies, explants and resections were each scored by the three pathologists. For ordinal variables, Light’s kappa (square weighted, for > 2 raters) was calculated using the 'psy'58 package, and intraclass correlation coefficient, Krippendorf’s alpha, and Kendall’s W were calculated using the 'irr'59 package. For the single scored categorical variable (histological diagnosis of NASH), Light’s kappa and Krippendorf’s alpha were calculated.
Each pathologist was randomly allocated a third of the cases and scored the whole slide images within the OMERO platform, entering the scores as key-value pairs associated with each image (NAS and SAF features for H&E, NASH-CRN and modified Ishak scores for PSR).
Digital pathology quantification of fat and fibrosis
A duplexed classification script was applied to the raw .ndpi scans of PSR-stained sections. An initial ‘whole tissue’ classifier was applied and small artefacts on the image were removed from the mask based on size. A second pixel-classifier of random trees (‘RTrees’) type with ‘gaussian’ and ‘weighted deviation’ features selected, pre-trained by one pathologist (TJK), was applied to classify pixels into the histological classes “fat”, “psr”, and “other tissue” within the “whole tissue” masked area, and the total number of pixels for each class and the percentage fat and PSR-positive pixels were returned. Whole slide scan processing was undertaken using command-line execution of a QuPath60 script on the University of Edinburgh’s Edinburgh Compute and Data Facility (ECDF) Linux compute cluster (Eddie). Classified images were reviewed by one pathologist (TJK), independent of all other data, and cases were excluded where large amounts of liver capsule or other artefact had been included in the classification.
Time-to-event analysis
Analysis using clinical and histopathological data only was undertaken in R61 (version 4.1.0) using the packages 'survival'62, 'survminer'63, and 'finalfit'64. Only needle biopsy cases were used for time-to-event analysis. Decompensation events in the clinical data extract were those defined by a combination of ICD codes and UK Operations/Procedure OPCS4 codes identifying activity relating to cirrhosis-related hospital admissions activity20,21. Decompensation and HCC event analysis was only undertaken on biopsy cases where the first decompensation or HCC-related coding was present in the clinical data extract after the biopsy date, and analysis was undertaken using death as a competing risk using Cox regression on the cause-specific hazard and Fine-Gray regression for proportional hazards modelling of the subdistribution hazard as a recommended comprehensive approach65. Kaplan-Meier estimator curves of all-cause mortality or event for assigned NASH-CRN fibrosis stages were compared by Log-Rank testing corrected for multiple comparisons using the Benjamini and Hochberg method in ‘survminer’62. Clustering of survival/event curves was determined using a k-means bootstrapped method19.Hazard ratios and 95% confidence intervals were derived from standard Cox regression models comparing mortality and rates of new-onset decompensation events with NASH-CRN stage, age at biopsy date, and gender as covariates. The assumption of proportional hazards for fitted Cox regression models was verified by examining plots of scaled Schoenfeld residuals against time for each covariate.
Multiplex immunofluorescence
Multiplexed immunofluorescence staining was performed using the MultiOmyx platform according to Gerdes et al66. This technology was performed using a single 4 µM FFPE slide where for each staining round two cyanine dye-labelled (Cy3, Cy5) antibodies were paired together. A custom multiplex panel was created consisting of 17 markers for a total of 9 rounds of antibody staining performed in sequence on the FFPE slides. In the case of TREM2 and CD9 which were applied as primary-secondary antibodies, samples were incubated with primary TREM2 and CD9 antibody followed by incubation with a species-specific secondary antibody conjugated to cyanine 5 or 3 (Cy5 or Cy3). The staining signal was then imaged and followed by a dye inactivation step, enabling repeated rounds of staining. The proprietary deep learning-based workflow NeoLYTX was subsequently applied to identify individual cells and perform cell classification for each marker and the phenotype of each cell was then determined through co-expression analysis.
Antibodies for MultiOmyx
Antibodies, by staining order, were rabbit anti-TREM2 (polyclonal, ProteinTech), mouse anti-MNDA (253A, Abcam), rabbit anti-CD9 (EPR2949, Abcam), mouse anti-CD66b (G10F5, BioLegend), mouse anti-CD11B (238439, R&D Systems), rabbit anti-DC-SIGN (D7F5C, Cell Signaling Technology), rabbit anti-Ki67 (SP6, Abcam), rabbit anti-IDO1 (SP260, Abcam), rabbit anti-CD11c (D3V1E, Cell Signaling), rabbit anti-PD-L1 (SP142, Abcam), rabbit anti-CD14 (EPR3652, Abcam), mouse anti-CD16 (DJ130c, ThermoFisher Scientific), mouse anti-CD68 (KP1, BioLegend), mouse anti-CD163 (EDHu-1, Bio-Rad), mouse anti-HLA DQ/DR/DP (WR18, Novus), mouse anti-CD33 (44M12D3, Novus Biologicals), mouse anti-SMA (1A4, Sigma-Aldrich).
RNA-sequencing methods
Hepatic RNA extraction and library preparation
RNA extraction was performed by the Genetics Core (Clinical Research Facility, Western General Hospital, Edinburgh, UK) from 4 x 10 µm curls of FFPE archival human NAFLD liver, using the QIAGEN miRNeasy FFPE Kit (Qiagen, UK) according to the manufacturer’s instructions. Quality control was by Qubit to measure RNA yield (and any potential DNA contamination) and Agilent Bioanalyser (Agilent, UK) to assess DV200. Samples with DV200 below 30% were not progressed for sequencing.
Libraries were prepared at PMS-IC using the low input Takara Bio SMARTer Stranded Total RNA-Seq Kit v2 (Pico Input Mammalian; Takara Bio, Sweden).
RNA sequencing and differential gene expression analysis
Sequencing data were generated by Edinburgh Genomics (University of Edinburgh, UK) using the Illumina NovaSeq 6000 platform (Illumina, Inc., US). Libraries were sequenced over a total of 22 S2 flow cells. Reads were trimmed using ‘Cutadapt’67 and aligned to the reference genome (GRCh38) using ‘STAR’68. Reads were assigned to features using 69 with a ’gtf’ file from Ensembl release 84.
For further analyses, R (version 4.1.2) was used. Samples with less than 1 million counts or 70% mapped reads were removed. Any specimen found to include HCC on histological review was excluded from further RNA-seq analyses.
Reads were normalized using the weighted trimmed mean of M-values (TMM) method70. Differential analysis was carried out with ‘limma-voom’ with the protein-coding genes. Statistical significance of genes was determined by an adjusted P-value according to Benjamini-Hochberg procedure of P < 0.05 and absolute FC of ≥ 171.
Principal component analysis (PCA) was performed to identify covariates that significantly correlated with the main principal components, so they could be controlled for downstream analyses (Extended Data Fig. 11). For this reason, sex was included as an additive effect in the linear model used for differential expression. Gene Set Enrichment Analysis (GSEA)72 was performed with GSEA function from ‘clusterProfiler’ (version 4.0.5)73. Data were visualized with ‘ggplot2’ (version 3.3.5)74 and ‘clusterProfiler’.
SNP variant calling
For genes PNPLA3, HSD17B13, TM6SF2 and MBOAT7 containing prespecified variants (“rs738409”, “rs72613567”, “rs58542926” and “rs641738” respectively), the Genome Analysis Toolkit (GATK)75 was used to call genotypes at all sites. The ‘HaplotypeCaller’ tool from GATK was used to produce a single genomic variant call format (GVCF) file per sample, with the following parameters: “–standard-min-confidence-threshold-for-calling 20”; “-dont-use-soft-clipped-bases"; “-L genes.bed”; “-ERC GVCF” where “genes.bed” is a bed file containing the genome coordinates of the genes of interest (padded by 1Kb in both directions to cover up- and down-stream SNPs). The ‘GenotypeGVCFs’ tool was then run with the parameter: “–include-non-variant-sites" to produce a single VCF file per sample. ‘SnpEff’ was used to functionally annotate the variants in the VCF files.
Deconvolution analysis
The MUlti-Subject SIngle Cell (‘MuSiC’) deconvolution tool45 was run using R version (3.6.3) to compare the counts from each bulk RNA-seq sample to reference single cell RNA-seq (scRNA-seq) data from five healthy and five cirrhotic livers46 to estimate the proportions of each sub-type of various lineages of non-parenchymal cells. This was done separately for each of the macrophage (non-circulating), B cells, mesenchyme, endothelial and T-cell lineages. In each analysis, only the scRNA-seq data from cells annotated as belonging to the lineage was used as a reference in the deconvolution, so the sum of estimated proportions of subtypes in each lineage equals one.
To assess the association between histological scores and proportions of various cell subtypes, the Spearman’s partial rank correlation coefficient was calculated using the R package ‘ppcor’ between the proportion of each subtype and the histological score, accounting for age and gender. Gender was converted to a numerical value by setting “male” to “1” and “female” to “2”. The histological data contained the “fibrosis stage”, which includes scores “1a”, “1b” and “1c”, with “1a” being the least severe, and “1c” being the most severe. To perform a rank correlation, these scores were converted to 1, 1.25 and 1.5 respectively. Association between proportion of each subtype and patient outcome (all-cause mortality or a decompensation event) was also assessed using Pearson’s correlation. For any given timepoint, samples were given a “0” in the event had not occurred, given a “2” if the event had occurred, and excluded if the patient was censored at that time.
Samples with less than one million reads assigned to genes in the bulk RNA-seq analysis, that lack a patient age, or do not come from a biopsy, were omitted.
Transcriptional risk score prediction model construction and performance assessment
DEGs between F3/F4 and F0/F1 biopsy cases were identified with ‘limma-voom’ in R. The criteria for statistical significance were adjusted P < 0.05 and absolute FC ≥ 1. Univariate Cox analysis was carried out to determine the DEGs associated with decompensation events according to an adjusted P < 0.01. Next, LASSO-penalised Cox regression was used to filter features using ‘glmnet’ package. To calculate a risk score, we used the following formula for each patient: \({\sum }_{i=1}^{n}{expression}_{i}* {coefficient}_{i}\). According to the median score, patients were divided into high- or low-risk groups. Kaplan-Meier analysis was performed using ‘survival’ and ‘survminer’ packages. A time-dependent receiver operating characteristic (ROC) curve was created by the ‘timeROC’ package to evaluate the predictive ability of the signature.
Transcriptional regulatory network construction and regulon analysis
Transcriptional network inference and regulon analysis was undertaken in R (version 4.1.0) using the ‘RTN’ package76. Normalised gene expression of all sequenced SteatoSITE samples was used to infer regulons corresponding to all TFs complied by Lambert et al30 by permutation analysis with the non-parametric estimator of mutual information and a P-value cut-off of 1.75e-7 to correct for multiple hypothesis testing. Unstable gene-TF interactions were removed by bootstrap analysis to produce a consensus network. The ARACNe algorithm77 using the data processing inequality theorem to enrich regulons with direct TF-gene target interactions was used to remove the weakest interaction in any triplet composed of two TFs with a common gene target; the threshold for filtering was determined from the null distribution derived in the permutation and bootstrapping steps.
To determine regulon activity scores for each sample, the tni.gsea2() function was used to calculate differential expression of each gene relative to the expression in the whole cohort and the ranked list of genes used to form a differential expression signature which was used alongside the TRN to determine differential sample-specific regulon activity.
Regulon activity profiles across disease stage were placed into 4 clusters by unsupervised soft clustering using the ‘Mfuzz’ package. The fuzzifier parameter, m, was directly estimated from the data using mestimate() to apply the method of Schwaemmle and Jensen78.
The ‘RTNsurvival’ package79 was used to undertake time-to-event analysis with THRB regulon activity. The cutpoint for THRB regulon activity was calculated using surv_cutpoint() from the ‘survminer’ package, applying maximally selected rank statistics of the ‘maxstat’ package with a minimum proportion of 0.25.
Shiny Interface
We developed an open-access gene browser, based on the START app80, that will ultimately be hosted on the SteatoSITE website (https://steatosite.com/). Shiny is an R package that allows the production of interactive plots and charts and is already widely used to allow users to explore RNA-seq data and create their own customised figures. A shiny interface was created to allow scientists to visualise the results of the analysis of differential gene expression according to NASH-CRN fibrosis stage.
Data availability
Hepatic bulk RNA-seq data is deposited in the European Nucleotide Archive (https://www.ebi.ac.uk/ena; study accession number: PRJEB58625). Gene expression data is also freely available for user-friendly interactive browsing online at https://shiny.igc.ed.ac.uk/SteatoSITE_gene_explorer/. SteatoSITE has delegated ethics allowing the granting of access to the full dataset within the PMS-IC secure environment to third parties by application (details at https://steatosite.com/researchers/), overseen and reviewed by the SteatoSITE Scientific Advisory Board.
Code availability
R scripts enabling the main steps of the analysis are available from the corresponding author on reasonable request.
Acknowledgements
We are grateful to the following for assistance: Lynn McMahon, Marian McNeil (PMS-IC) - project management, data compliance and security, workspace provision; Marc Jones (PMS-IC) - RNA-seq libraries; Emily Coulston (PMS-IC) - administrative support; Edinburgh Clinical Research Facility Genetics Core/Lee Murphy – RNA extraction; Edinburgh Genomics/Matt Arno – RNA-seq; electronic Data Research and Innovation Service (eDRIS)/Beth Bruce – clinical data access approval; Safe Havens: Lothian (Chris Duncan, Pamela Linksted), Grampian (Katie Wilde, Joanne Lumsden), GG&C (Jennie Burnie, Samuel McKenna) - EHR data retrieval; Lothian Bioresource (Linda McLeod, Craig Marshall) - human liver tissue access; Ewan McDowall (Institute of Genetics and Cancer) – Shiny app hosting.
This study was funded by Innovate UK (Precision medicine: impacting through innovative technology (Reference: TS/R017581/1)). MJR is an iCASE PhD student funded by the Medical Research Council Precision Medicine Doctoral Training Programme (Reference: MR/R01566X/1) and Galecto Biotech. FM is supported by an EPSRC Innovation Fellowship (Reference: EP/S001921/1), an EPSRC New Investigator Award (Reference: EP/R035350/1), and a UKRI Engineering Biology Transition Award (Reference: BB/W014610/1). LB is supported by Wellcome-University of Edinburgh Institutional Strategic Support Fund (Reference: 204804/Z/16/Z) and an EPSRC Postdoctoral Fellowship (Reference: EP/P017134/1).
Author information
Affiliations
Institute for Regeneration and Repair, University of Edinburgh, Edinburgh, UK
T.J. Kendall, M. J. Ramos, P. Ramachandran, J. A. Fallowfield
Edinburgh Pathology, University of Edinburgh, Edinburgh, UK
T.J. Kendall
Precision Medicine Scotland-Innovation Centre (PMS-IC), University of Glasgow, Glasgow, UK
M. McColgan, H. Ellis
PMS-IC is a consortium consisting of four NHS Scotland health boards (Lothian, Tayside, Grampian, and Greater Glasgow & Clyde) and the four medical universities of Scotland (Edinburgh, Glasgow, Aberdeen, and Dundee), in addition to two companies (Aridhia and Thermo Fisher).
Edinburgh Genomics (Bioinformatics), University of Edinburgh, Edinburgh, UK
F. Turner, D. Dunbar
OHSU-PSU School of Public Health, Oregon Health & Sciences University, Portland, OR 97239, USA
J. Minnier
Knight Cancer Institute Biostatistics Shared Resource, Oregon Health & Sciences University, Portland, OR 97239, USA
J. Minnier
Pathology Department, Queen Elizabeth University Hospital, Glasgow, UK
G. Kohnen, P. Konanahalli, K. Oien
National Institute of Health Research (NIHR) Nottingham Biomedical Research Centre, Nottingham University Hospitals NHS Trust and University of Nottingham, Nottingham, UK
I.N. Guha
School of Engineering, Institute of Bioengineering, University of Edinburgh, Edinburgh, UK
L. Bandiera, F. Menolascina
Centre for Engineering Biology, University of Edinburgh, Edinburgh, UK
L. Bandiera, F. Menolascina
NeoGenomics Laboratories, Fort Myers, Florida, USA
A. Juncker-Jensen
NHS Greater Glasgow and Clyde Safe Haven, Glasgow, UK
C. Mayor, D. Alexander
Contributions
JAF, TJK and MDM contributed to the development of the study concept and design and JAF, TJK, MJR, FT, HE, GK, PK, PR, MDM, DD, LB, FM, AJ-J, JM and ING contributed to data analysis and interpretation. JAF and TJK drafted the initial manuscript. TJK, GK, PK, KO, CM, DA, DD and MDM provided technical and/or material support for the project. JAF is the guarantor of the work. Each author contributed important intellectual content during manuscript drafting or revision and accepts accountability for the overall work by ensuring that questions pertaining to the accuracy or integrity of any portion of the work are appropriately investigated and resolved. All authors approved the final version of the manuscript. The corresponding author attests that all the listed authors meet the authorship criteria and that no others meeting the criteria have been omitted.
Corresponding author
Correspondence to J. A. Fallowfield ([email protected]).
Competing interests
TJK serves as a consultant for, or has received speakers' fees from, Resolution Therapeutics, Clinovate Health, Perspectum Diagnostics, and Incyte Corporation. PR serves as a consultant for Merck and has received research grant funding from Genentech. JAF serves as a consultant or advisory board member for Ipsen, Redx Pharma, River 2 Renal Corp., Stimuliver, Galecto Biotech, ICON and Guidepoint, and has received research grant funding from Intercept Pharmaceuticals and Genentech. AJJ is an employee and stock owner at NeoGenomics. ING serves as an advisory board member for Resolution Therapeutics and has received research grant funding from Gilead Sciences. MDM and DD have a controlling shareholder interest in Biodev Ltd.
- WHO Regional Office for Europe. ‘SDR, chronic liver disease and cirrhosis, all ages, per 100 000’ European Health for All explorer. https://gateway.euro.who.int/en/indicators/hfa_236-1860-sdr-chronic-liver-disease-and-cirrhosis-all-ages-per-100-000/ (2021).
- Younossi, Z. M. et al. The economic and clinical burden of nonalcoholic fatty liver disease in the United States and Europe. Hepatology 64, 1577–1586 (2016).
- NHS Blood and Transplant. NHS Blood and Transplant Annual Report and Accounts 2018/19. (2019).
- Le, M. H. et al. Forecasted 2040 global prevalence of nonalcoholic fatty liver disease using hierarchical bayesian approach. Clin Mol Hepatol 28, 841–850 (2022).
- Loomba, R., Friedman, S. L. & Shulman, G. I. Mechanisms and disease consequences of nonalcoholic fatty liver disease. Cell 184, 2537–2564 (2021).
- Cai, J., Zhang, X.-J. & Li, H. The Role of Innate Immune Cells in Nonalcoholic Steatohepatitis. Hepatology 70, 1026–1037 (2019).
- Anstee, Q. M. et al. Genome-wide association study of non-alcoholic fatty liver and steatohepatitis in a histologically characterised cohort☆. J Hepatol 73, 505–515 (2020).
- Loomba, R. et al. Gut Microbiome-Based Metagenomic Signature for Non-invasive Detection of Advanced Fibrosis in Human Nonalcoholic Fatty Liver Disease. Cell Metab 25, 1054-1062.e5 (2017).
- Meijnikman, A. S. et al. Microbiome-derived ethanol in nonalcoholic fatty liver disease. Nat Med 28, 2100–2106 (2022).
- Angulo, P. et al. Liver Fibrosis, but no Other Histologic Features, Associates with Long-term Outcomes of Patients With Nonalcoholic Fatty Liver Disease. Gastroenterology 149, 389-397.e10 (2015).
- Vilar-Gomez, E. et al. Fibrosis Severity as a Determinant of Cause-Specific Mortality in Patients With Advanced Nonalcoholic Fatty Liver Disease: A Multi-National Cohort Study. Gastroenterology 155, 443-457.e17 (2018).
- Sanyal, A. J. et al. Prospective Study of Outcomes in Adults with Nonalcoholic Fatty Liver Disease. New England Journal of Medicine 385, 1559–1569 (2021).
- Rowe, I. A. & Parker, R. The Placebo Response in Randomized Trials in Nonalcoholic Steatohepatitis Simply Explained. Clin Gastroenterol Hepatol 20, e564–e572 (2022).
- Ratziu, V. & Friedman, S. L. Why do so many NASH trials fail? Gastroenterology S0016-5085(20)30680–6 (2020) doi:10.1053/j.gastro.2020.05.046.
- Chopra, H., Baig, A. A., Gautam, R. K. & Kamal, M. A. Application of Artificial intelligence in Drug Discovery. Curr Pharm Des (2022) doi:10.2174/1381612828666220608141049.
- Asiimwe, R. et al. From biobank and data silos into a data commons: convergence to support translational medicine. J Transl Med 19, 493 (2021).
- Kleiner, D. E. et al. Design and validation of a histological scoring system for nonalcoholic fatty liver disease. Hepatology 41, 1313–1321 (2005).
- Jung, E. S. et al. Interobserver Agreement on Pathologic Features of Liver Biopsy Tissue in Patients with Nonalcoholic Fatty Liver Disease. JPTM 50, 190–196 (2016).
- Villanueva, N. M., Sestelo, M. & Meira-Machado, L. A method for determining groups in multiple survival curves. Statistics in Medicine 38, 866–877 (2019).
- Hagström, H. et al. Administrative Coding in Electronic Health Care Record-Based Research of NAFLD: An Expert Panel Consensus Statement. Hepatology 74, 474–482 (2021).
- Innes, H. et al. Performance of routine risk scores for predicting cirrhosis-related morbidity in the community. Journal of Hepatology 77, 365–376 (2022).
- Piscaglia, F. et al. Clinical patterns of hepatocellular carcinoma in nonalcoholic fatty liver disease: A multicenter prospective study. Hepatology 63, 827–838 (2016).
- Stine, J. G. et al. Systematic review with meta-analysis: risk of hepatocellular carcinoma in non-alcoholic steatohepatitis without cirrhosis compared to other liver diseases. Alimentary Pharmacology & Therapeutics 48, 696–703 (2018).
- Illumina. Evaluating RNA Quality from FFPE Samples. (2015).
- Govaere, O. et al. Transcriptomic profiling across the nonalcoholic fatty liver disease spectrum reveals gene signatures for steatohepatitis and fibrosis. Science Translational Medicine 12, eaba4448 (2020).
- Suppli, M. P. et al. Hepatic transcriptome signatures in patients with varying degrees of nonalcoholic fatty liver disease compared with healthy normal-weight individuals. American Journal of Physiology-Gastrointestinal and Liver Physiology 316, G462–G472 (2019).
- Ramachandran, P. et al. Resolving the fibrotic niche of human liver cirrhosis at single cell level. Nature 575, 512–518 (2019).
- Hu, X. et al. Polygenic transcriptome risk scores for COPD and lung function improve cross-ethnic portability of prediction in the NHLBI TOPMed program. American Journal of Human Genetics 109, 857 (2022).
- Liang, Y. et al. Polygenic transcriptome risk scores (PTRS) can improve portability of polygenic risk scores across ancestries. Genome Biol 23, 23 (2022).
- Lambert, S. A. et al. The Human Transcription Factors. Cell 172, 650–665 (2018).
- Carro, M. S. et al. The transcriptional network for mesenchymal transformation of brain tumours. Nature 463, 318–325 (2010).
- Ng, C. H. et al. Mortality Outcomes by Fibrosis Stage in Nonalcoholic Fatty Liver Disease: A Systematic Review and Meta-analysis. Clinical Gastroenterology and Hepatology 0, (2022).
- Ramachandran, P., Matchett, K. P., Dobie, R., Wilson-Kanamori, J. R. & Henderson, N. C. Single-cell technologies in hepatology: new insights into liver biology and disease pathogenesis. Nat Rev Gastroenterol Hepatol 17, 457–472 (2020).
- Chu, A. L., Schilling, J. D., King, K. R. & Feldstein, A. E. The Power of Single Cell Analysis for the Study of Liver Pathobiology. Hepatology 73, 437–448 (2021).
- Wallace, S. J., Tacke, F., Schwabe, R. F. & Henderson, N. C. Understanding the cellular interactome of non-alcoholic fatty liver disease. JHEPReport 4, (2022).
- Buisseret, L. et al. Clinical significance of CD73 in triple-negative breast cancer: multiplex analysis of a phase III clinical trial. Ann Oncol 29, 1056–1062 (2018).
- Bosisio, F. M. et al. Next-Generation Pathology Using Multiplexed Immunohistochemistry: Mapping Tissue Architecture at Single-Cell Level. Front Oncol 12, 918900 (2022).
- Ilie, M. et al. Automated chromogenic multiplexed immunohistochemistry assay for diagnosis and predictive biomarker testing in non-small cell lung cancer. Lung Cancer 124, 90–94 (2018).
- D’Amico, G., Garcia-Tsao, G. & Pagliaro, L. Natural history and prognostic indicators of survival in cirrhosis: a systematic review of 118 studies. J Hepatol 44, 217–231 (2006).
- Bianco, C. et al. Non-invasive stratification of hepatocellular carcinoma risk in non-alcoholic fatty liver using polygenic risk scores. J Hepatol 74, 775–782 (2021).
- Vincentis, A. D. et al. A Polygenic Risk Score to Refine Risk Stratification and Prediction for Severe Liver Disease by Clinical Fibrosis Scores. Clinical Gastroenterology and Hepatology 20, 658–673 (2022).
- Marigorta, U. M. et al. Transcriptional Risk Scores link GWAS to eQTL and Predict Complications in Crohn’s Disease. Nat Genet 49, 1517–1521 (2017).
- Wu, I. & Moses, M. A. BNF-1, a novel gene encoding a putative extracellular matrix protein, is overexpressed in tumor tissues. Gene 311, 105–110 (2003).
- Huang, G. & Brigstock, D. R. Regulation of hepatic stellate cells by connective tissue growth factor. Front Biosci (Landmark Ed) 17, 2495–2507 (2012).
- Chan, K. K.-S. et al. Stanniocalcin 1 is a serum biomarker and potential therapeutic target for HBV-associated liver fibrosis. J Pathol 257, 227–238 (2022).
- Tao, L. et al. Glial cell line-derived neurotrophic factor (GDNF) mediates hepatic stellate cell activation via ALK5/Smad signalling. Gut 68, 2214–2227 (2019).
- Liu, Z. et al. Wnt-TCF7-SOX9 axis promotes cholangiocarcinoma proliferation and pemigatinib resistance in a FGF7-FGFR2 autocrine pathway. Oncogene 41, 2885–2896 (2022).
- Wirth, E. K., Puengel, T., Spranger, J. & Tacke, F. Thyroid hormones as a disease modifier and therapeutic target in nonalcoholic steatohepatitis. Expert Rev Endocrinol Metab 17, 425–434 (2022).
- Harrison, S. A. et al. Effects of Resmetirom on Noninvasive Endpoints in a 36-Week Phase 2 Active Treatment Extension Study in Patients With NASH. Hepatol Commun 5, 573–588 (2021).
- The Scottish Public Health Observatory. Chronic liver disease: international comparisons. https://www.scotpho.org.uk/health-wellbeing-and-disease/chronic-liver-disease/data/international-comparisons/ (2022).
- Shearer, J. E. et al. Systematic review: development of a consensus code set to identify cirrhosis in electronic health records. Aliment Pharmacol Ther 55, 645–657 (2022).
- Staufer, K. et al. Ethyl glucuronide in hair detects a high rate of harmful alcohol consumption in presumed non-alcoholic fatty liver disease. J Hepatol 77, 918–930 (2022).
- Sookoian, S., Castaño, G. O. & Pirola, C. J. Modest alcohol consumption decreases the risk of non-alcoholic fatty liver disease: a meta-analysis of 43 175 individuals. Gut 63, 530–532 (2014).
- Jarvis, H. et al. Does moderate alcohol consumption accelerate the progression of liver disease in NAFLD? A systematic review and narrative synthesis. BMJ Open 12, e049767 (2022).
- Allan, C. et al. OMERO: flexible, model-driven data management for experimental biology. Nat Methods 9, 245–253 (2012).
- Bedossa, P. et al. Histopathological algorithm and scoring system for evaluation of liver lesions in morbidly obese patients. Hepatology 56, 1751–1759 (2012).
- McDonald, N. et al. Multiparametric magnetic resonance imaging for quantitation of liver disease: a two-centre cross-sectional observational study. Sci Rep 8, (2018).
- Falissard, B. psy: Various Procedures Used in Psychometrics. (2022).
- Gamer, M., Lemon, J. & Singh <[email protected]>, I. F. P. irr: Various Coefficients of Interrater Reliability and Agreement. (2019).
- Bankhead, P. et al. QuPath: Open source software for digital pathology image analysis. Scientific Reports 7, 16878 (2017).
- R Core Team. R: A Language and Environment for Statistical Computing. (R Foundation for Statistical Computing, 2022).
- Kassambara, A., Kosinski, M. & Biecek, P. survminer: Drawing Survival Curves using ‘ggplot2’. (2021).
- Therneau, T. M. A Package for Survival Analysis in R. (2022).
- Harrison, E., Drake, T. & Ots, R. finalfit: Quickly Create Elegant Regression Results Tables and Plots when Modelling. (2022).
- Latouche, A., Allignol, A., Beyersmann, J., Labopin, M. & Fine, J. P. A competing risks analysis should report results on all cause-specific hazards and cumulative incidence functions. J Clin Epidemiol 66, 648–653 (2013).
- Gerdes, M. J. et al. Highly multiplexed single-cell analysis of formalin-fixed, paraffin-embedded cancer tissue. Proc Natl Acad Sci U S A 110, 11982–11987 (2013).
- Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.journal 17, 10–12 (2011).
- Dobin, A. et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21 (2013).
- Liao, Y., Smyth, G. K. & Shi, W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30, 923–930 (2014).
- Robinson, M. D. & Oshlack, A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biology 11, R25 (2010).
- Storey, J. D. & Tibshirani, R. Statistical significance for genomewide studies. Proceedings of the National Academy of Sciences 100, 9440–9445 (2003).
- Subramanian, A. et al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proceedings of the National Academy of Sciences 102, 15545–15550 (2005).
- Wu, T. et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb) 2, 100141 (2021).
- Wickham, H. ggplot2: Elegant Graphics for Data Analysis. (Springer-Verlag New York, 2016).
- Van der Auwera, G. A. et al. From FastQ data to high confidence variant calls: the Genome Analysis Toolkit best practices pipeline. Curr Protoc Bioinformatics 43, 11.10.1-11.10.33 (2013).
- Castro, M. A. A. et al. Regulators of genetic risk of breast cancer identified by integrative network analysis. Nat Genet 48, 12–21 (2016).
- Margolin, A. A. et al. Reverse engineering cellular networks. Nat Protoc 1, 662–671 (2006).
- Schwämmle, V. & Jensen, O. N. A simple and fast method to determine the parameters for fuzzy c–means cluster analysis. Bioinformatics 26, 2841–2848 (2010).
- Groeneveld, C. S. et al. RTNsurvival: an R/Bioconductor package for regulatory network survival analysis. Bioinformatics 35, 4488–4489 (2019).
- Nelson, J. W., Sklenar, J., Barnes, A. P. & Minnier, J. The START App: a web-based RNAseq analysis and visualization resource. Bioinformatics 33, 447–449 (2017).
Table 1. Demographics and clinical characteristics of SteatoSITE cohort.
Yes there is potential Competing Interest.
TJK serves as a consultant for, or has received speakers' fees from, Resolution Therapeutics, Clinovate Health, Perspectum Diagnostics, and Incyte Corporation. PR serves as a consultant for Merck and has received research grant funding from Genentech. JAF serves as a consultant or advisory board member for Ipsen, Redx Pharma, River 2 Renal Corp., Stimuliver, Galecto Biotech, ICON and Guidepoint, and has received research grant funding from Intercept Pharmaceuticals and Genentech. AJJ is an employee and stock owner at NeoGenomics. ING serves as an advisory board member for Resolution Therapeutics and has received research grant funding from Gilead Sciences. MDM and DRD have a controlling shareholder interest in Biodev Ltd.
- SupplementaryTables.xlsx
- editorialpolicychecklistcomplete.pdf
Editorial Policy Checklist
- reportingsummarycomplete.pdf
Reporting Summary
- SupplementaryInformation.docx