A robust gene expression signature for NASH in liver expression data

Non-Alcoholic Fatty Liver Disease (NAFLD) is a progressive liver disease that affects up to 30% of worldwide population, of which up to 25% progress to Non-Alcoholic SteatoHepatitis (NASH), a severe form of the disease that involves inflammation and predisposes the patient to liver cirrhosis. Despite its epidemic proportions, there is no reliable diagnostics that generalizes to global patient population for distinguishing NASH from NAFLD. We performed a comprehensive multicohort analysis of publicly available transcriptome data of liver biopsies from Healthy Controls (HC), NAFLD and NASH patients. Altogether we analyzed 812 samples from 12 different datasets across 7 countries, encompassing real world patient heterogeneity. We used 7 datasets for discovery and 5 datasets were held-out for independent validation. Altogether we identified 130 genes significantly differentially expressed in NASH versus a mixed group of NAFLD and HC. We show that our signature is not driven by one particular group (NAFLD or HC) and reflects true biological signal. Using a forward search we were able to downselect to a parsimonious set of 19 mRNA signature with mean AUROC of 0.98 in discovery and 0.79 in independent validation. Methods for consistent diagnosis of NASH relative to NAFLD are urgently needed. We showed that gene expression data combined with advanced statistical methodology holds the potential to serve basis for development of such diagnostic tests for the unmet clinical need.

In this study, we hypothesized that a multi-cohort analysis of transcriptome profiles of liver biopsies from patients with NAFLD or NASH would identify a robust signature for NASH that generalizes across the biological, clinical, and technical heterogeneity of the real-world patient population and will be suitable for clinical development.

Results
Data collection and compilation. Our search in GEO and ArrayExpress for transcriptome profiles of liver biopsies from patients with NASH or NAFLD resulted in 12 datasets composed of 812 samples from patients across 7 countries that met our inclusion criteria (Table 1, methods). Overall, the 12 datasets represent a broad spectrum of biological, clinical, and technical heterogeneity. They include samples from adolescents or adults, with or without comorbidities, and profiled using a diverse set of commercial microarrays and RNA sequencing platforms. HCs across these datasets represented real-world heterogeneity as they included those with normal weight, healthy obese, and those suspected of NAFLD. We used 309 samples from 7 datasets identifying differentially expressed genes and 503 samples from 5 datasets for independent validation of those genes. When selecting datasets for discovery and validation, we aimed to maximize the biological, clinical, and technical heterogeneity representation within the discovery datasets, while ensuring that the number of samples allocated to discovery to no more than 50% of the total samples. We first identified differentially expressed genes in NASH compared to NAFLD and HC groups as our primary comparison. We subsequently extended our analysis to all 6 possible comparisons (Fig. 1A, Table 2). Notably, not all datasets included samples from all phenotypic groups, thus some datasets were excluded from particular comparisons. Membership of datasets in particular comparisons is summarized in Table S1.
Gene expression signature differentiates NASH from NAFLD or healthy controls. First, we performed a multicohort analysis of NASH vs "others" where others included both the HC and NAFLD groups ([NASH]vs[NAFLD + HC]). At thresholds of |ES|> 0.6 and FDR < 10%, we identified 130 genes that were consistently differentially expressed in NASH vs others (85 over-and 45 under-expressed, Table 2). MetaIntegrator score derived from these 130 genes yielded mean discovery AUROC of 0.94 and validation AUROC of 0.80 ( Fig. 1B,C,D), and showed low inter-study heterogeneity of ES for particular genes (e.g. FABP4, Fig. 1E). This strongly suggests that our signature is likely to generalize well in independent datasets, and that a true and robust gene expression signal can be harnessed to reliably distinguish NASH from less severe forms of NAFLD in real world patient populations.
Despite the encouraging performance of the [NASH]vs[NAFLD + HC] signature, it is possible that widespread study design or sample ascertainment methods may have introduced technical variance between the NASH, NAFLD and HC groups, thus artificially inflating any signal in our data. To evaluate this possibility, we explored the performance of all other 5 possible signatures in this data (Fig. 1A,E). For a true biological signal we would expect: (a) signature [NASH]vs[HC] to carry the clearest signal and perform the best; (b) the signature [NAFLD]vs[HC] to carry a weaker signal and perform worse than any signature that involves NASH; and (c) given the progressive nature of the disease [NAFLD]vs[NASH + HC] is unlikely to produce a meaningful signal, serving as an "internal scramble control". Conversely, if there is a persistent artificial difference between the groups introduced by the sampling process, we would expect that whatever comparison we make would validate similarly well. Indeed, the results from the other 5 signatures (Fig. 1E)  Gene expression differences validate known pathways and suggests novel biology. To gain biological insight we compared the differentially expressed genes between the signatures. Each signature has between 41 and 173 differentially expressed genes ( Table 2) and altogether the six signatures encompass 428 genes (Table S1). The four signatures that include NASH were consistent with each other. For example, effect sizes (ESs) for genes captured in those four contrasts maintain their directionality ( Fig. 2A) and the four signatures that include NASH patients as cases share a substantial number of genes among them ( Fig. 2A     www.nature.com/scientificreports/ from NAFLD or HCs, measuring expression of 130 genes in clinical setting might be unpractical and unnecessary. In addition, while all the genes we used are differentially expressed in the NASH group, considerable correlation between the expression of some signature genes renders their contribution to the performance of the final score redundant. We used a greedy forward search to narrow down the number of genes necessary for our signature performance. We were able to minimize our gene set to 19 genes (Table S1) while optimizing for discriminatory performance between NASH and NAFLD or HC groups in the 7 studies used for discovery. We then re-computed a diagnostic score using only those 19 genes, and show that the 19 gene score performs as well as the complete list (validation summary AUROC 0.79 vs summary AUROC 0.80; Fig. 3A,B).

Discussion
To identify a persistent and reproducible gene expression signature for NASH, we performed the most comprehensive meta-analysis of published liver gene expression to date, encompassing 12 datasets (812 samples) that represent real world patient populations. We successfully leveraged both biological and technical heterogeneity in identifying a gene expression signature that robustly distinguishes NASH from NAFLD or HCs. Multiple tests and scoring systems have been developed for non-invasive diagnostics for NAFLD and NASH. These perform well enough for detecting NAFLD, but not NASH, or are limited to particular populations 12,16 . Our 19-gene signature achieved a performance of mean AUROC = 0.98 in discovery and AUROC = 0.79 in independent, blind validation in 5 datasets that included multiple populations or phenotypes. This strongly suggests that our approach has great potential for development of gene-expression-based diagnostic test for NASH. Gene enrichment analysis of the NASH signature reflects the diagnostic hallmark of NASH, relative to NAFLD-inflammation. Interestingly, many of the genes we identify were previously implicated as potential markers for HCC survival and progression. This result suggests that the molecular processes involved in progression of NASH to cirrhosis and HCC are evident in at least some NASH patients and could serve as basis for development of endotyping strategy for NASH.
As part of this study, we also examined all six possible signatures, and show that the discriminatory performance of [ Ultimately, not all six signatures perform equally well in our analysis. It is clear that the distinction between NAFLD and HC ([NAFLD]vs[HC]) proved to be more difficult than distinguishing NASH from any combination of HC and NAFLD. This is expected as it mirrors the situation in the clinic. Also, the large spectrum of HC  (Table 1) suggests a significant overlap between definitions of NAFLD, healthy obese, and healthy. For example, HC samples in GSE66676 were taken from obese patients undergoing bariatric surgery. We believe that this wide variation in control population, at least in part, explains the relatively large variation in our validation performance. In particular, GSE83452, the worst performing validation dataset, is the only validation study that does not include patients diagnosed as NAFLD and uses obese as control population. The control samples in this dataset were taken from people who had abnormal liver enzymes and had a biopsy taken, which was declared not NAFLD by a pathologist. Thus, in our opinion, the relatively large drop in performance between the discovery and validation is at least in part due to the differences in diagnostic approach and potentially inter-individual differences in histological evaluations between different pathologists 18,55 . This variability further underscores the need for more cohesive and objective diagnostic framework for NASH. There are some limitations to our study. First, while our approach accounts for biological, technical, and clinical heterogeneity, it is limited by the variability captured in the datasets available to us. It is possible that incorporation of additional sources of variability-such as age, ethnicities, geographical areas, or technical platforms, would lead to reduced performance in these settings. Thus, further retrospective validation in independent cohorts is needed to ensure that we can refine our gene signature and achieve a level of performance that has clinical utility. Yet this also highlights the strength of our framework: as new datasets become available, they can and should be incorporated in the analysis producing updated and refined signature models. Second, to reach clinical utility, our 19-gene signature should undergo further algorithmic refinement. Notwithstanding, the computational refinement process can also be computationally expensive, presents unique challenges and ideally should be done when more data is available. Thus, we believe that this type of signature optimization is beyond the scope of this manuscript. Third, the holy grail of NASH diagnostic is non-invasive molecular test that would not require liver biopsy. Our analysis indicates that the [NASH]vs[NAFLD + HC] signature is driven by changes in composition of immune cells (e.g. leukocyte migration), and we expect to be able to detect such processes in blood as well. Our experience with other diseases 28,33,40,41,[56][57][58] suggests that given enough data, the MetaIntegrator framework would be successful in developing blood based signature. However, transcriptomic blood data for NASH patients is very sparce. Hence, our strategy is first to identify and validate NASH gene signature based on liver data and translate it into blood, when blood transcriptomic data becomes available. Finally, in this work we focused on feature selection which is only a first step in building a model. Our score is a difference of geometric means, it does not produce probability score and therefore some calibration metrics, such as calibration curves or Hosmer-Lemeshow test are not applicable. Rather, these metrics if applicable, should be included in final model evaluation.
Taken together, we believe that our work provides a solid foundation for development of gene-expressionbased test for NASH. Accurate liver gene expression testing could help with NASH diagnosis, patient care, and potentially drug development. For example, it could inform retrospective analysis of clinical trials of failed NASH therapies-potentially opening an avenue for more successful characterization of patient subgroups that benefited from the treatment and repositioning of otherwise failed drug candidates. We envision that, in addition to the current standards of histological evaluation by qualified pathologists, a gene-expression-based NASH diagnostic will add value in clinical decision making and promote standardization in the field. Biopsies represent only a small area of the liver, therefore histological changes that have not widely spread could be missed. In principle, a gene-expression-based diagnosis would have the added value of reflecting changes in cellular microenvironment that occur outside of that particular area.
To summarize, our results demonstrate that gene expression analysis harbors an exciting opportunity for development of diagnostic test for NASH. We believe that this work provides a solid foundation for further development, both in terms of algorithmic refinement of the presented signatures and addition of other datasets that would help develop accessible, high throughput and reliable diagnostic. With further prospective validation, our results hold the potential for breakthrough diagnostic test for NASH.

Methods
Data collection. We searched public gene expression repositories (GEO 59 and ArrayExpress 60 ) for datasets that included transcriptome profiles of liver biopsies from patients with NASH in January 2020. We identified 83 datasets with 3,359 samples. We excluded datasets that did not meet the following criteria: human, liver tissue, includes at least 5 patients in either NAFLD or HC groups and at least 5 patients in the NASH group. Twelve datasets met the inclusion criteria (Table 1). We then manually curated the 12 datasets to ensure integrity of phenotypic data, diagnostic criteria of NASH, NAFLD, and HC patients, and for general match between the deposited data and the numbers cited in associated publication when available.
Data preprocessing. For each dataset, we downloaded raw expression data and pre-processed using standard methods. Specifically, we applied RMA to all data from Affymetrix platforms 61,62 and used limma package 63,64 with quantile normalization for Illumina and other commercial arrays. For RNAseq datasets, we downloaded the associated SRA read fastq files and used FastQC for initial quality control. We then used STAR v2.2 65,66 , human genome GRCh38 67 and GENCODE 68,69 v32 human genome annotation for read alignment and gene expression quantification, as previously described 70 . To facilitate integrated analysis, we used the Annotation Dbi and Hs.org 71 packages to map probe and gene identifiers in each dataset to Entrez Gene identifiers (IDs). We used sample identifier to match phenotypic data from the databases to expression data and used the phenotypic and expression data to ensure sample uniqueness. We found no duplicated samples between the datasets. Multicohort analysis. We used the R package MetaIntegrator for multi-cohort analysis 26,27,72  We a priori divided the datasets into two groups: (1) 7 datasets (309 samples) as "discovery cohorts" and (2) 5 datasets (503 samples) held out as independent "validation cohort" (Table 1). When dividing the studies into discovery and validation we sought to maximize the technical and biological heterogeneity encompassed by the discovery studies, while keeping the number of samples for discovery bellow 50%. To avoid overfitting we also designate all studies from the same research group either for discovery or validation. For each comparison, we defined the more severe phenotypic group as "case" and the less severe group as "control" ( Table 2). We used the more stringent Leave-One-Study-Out (LOSO) cross-validation within the discovery cohorts to identify differentially expressed genes, whereby each of the discovery datasets is left out in a round-robin fashion to obtain ES for each gene. In each iteration, we applied the ES and q-value thresholds, and selected the genes that met these thresholds in every iteration. This ensures that gene's ES is not driven by one particular dataset and results in a list of genes that are consistently positively or negatively differentially expressed between the classes. We examined the gene lists produced by MetaIntegrator filterGenes function over multiple ES and q-value thresholds (0.6-0.8, and 0.1-0.01 respectively), adopting the cutoffs of |ES|> = 0.6 and FDR p < 0.1. Notably, the FDR in this case refers to the significance level associated with the pooled effect size, not to a separate testing of differential expression. These thresholds were determined based on the guideline from power estimate ( Figure S1) as well as to allow for meaningful pathway analysis in all 6 signatures. Expression of the selected genes is then combined into a score: Gene set enrichment analysis. We used enrichGO function from the R package clusterProfiler to perform gene enrichment analysis [73][74][75] . The package supports overrepresentation test against the entirety of organism specific GO annotation as represented in the OrgDb object, and provides Benjamini-Hochberg adjusted p-value for the observed overrepresentations. We used union of all genes expressed in the discovery datasets, for general background.
Approval for human experiments. All studies included in this meta-analysis obtained informed consent of the human subjects and were performed in accordance with relevant named guidelines and regulations governing the study. No new subjects were recruited for the purpose of this work.