1. Study search and selection

This review was conducted in June 2020, following the Preferred Reporting Items for Systematic Reviews and Meta-Analyses (PRISMA) statement guidelines (1). We searched the GEO datasets database (2) using the keywords: “NAFLD”, “NAFL”, “steatosis”, “NASH”, and “steatohepatitis” for transcriptomic studies published in English.

Study exclusion criteria

We applied the following exclusion criteria: (i) studies conducted in organisms other than humans; (ii) studies without information about the sex of the participants or that had not included both sexes; (iii) studies without individuals from the NAFL and NASH stages of the disease; (iv) studies in which the disease had not been diagnosed with a biopsy. In the latter case, we decided to require the use of a biopsy as the method for NAFLD and NASH diagnosis because other less invasive methods such as conventional imaging techniques have a low sensitivity to detect mild steatosis and/or to differentiate NAFL from NASH (3).

2. Bioinformatics analysis strategy

We applied the same three-step transcriptomic analysis strategy to each separate study: (i) data acquisition and preprocessing; (ii) differential gene expression and functional enrichment analysis; (iii) differential pathway activation and functional analysis. The functional results of all the studies were then integrated using meta-analysis techniques (Figure 1). This work was carried out using the programming language R 3.6.0 (4); information about the packages we used and their version numbers are provided in the supplementary material.

3. Data acquisition and preprocessing

During the data preprocessing we standardised the nomenclature of the different stages of the disease (by grouping the patients into Control, NAFL, and NASH groups) as well as the association of the probe identifiers with their corresponding genes. In those studies, with treated and untreated patients, treated patients were eliminated for further analysis. For studies that had used microarrays, the probe codes were transformed to their respective Entrez identifiers from the NCBI database (5).

For repeated probes, we calculated the median of their expression values. When available, data normalised by the original authors of the studies were used. Otherwise (studies GSE61260 and GSE66676), the raw data was normalised using the Robust Multichip Average (RMA) algorithm (6). The gene nomenclature in RNA-Seq studies was also standardised to the Entrez identifiers and the raw data matrix was processed using the TMM standardisation method (7) followed by a log-transformation of the data.

To meet the objective of this work, women were separated into premenopausal and postmenopausal groups. Since the selected studies did not include this information, we assumed that women aged under 50 years had been premenopausal, based on previous literature indicating that the average age of menopause is around 48–52 years (8,9).

After the normalisation of the data we attempted to detect any possible anomalous effects within the studies by completing an exploratory clustering and a principal component analysis (PCA).

4. Differential gene expression and functional enrichment analysis

We carried out a differential expression analysis of these studies using the limma and edgeR packages (10,11) to detect genes that were differentially expressed during the progression of NAFL to NASH in premenopausal women versus men. For each gene, we adjusted a linear model, which included possible batch effects, contrasting (NASH.W – NAFL.W) and (NASH.M – NAFL.M), where NASH.W, NAFL.W, NASH.M, and NAFL.M corresponded to NASH-affected premenopausal women, NAFL-affected premenopausal women, NASH-affected men, and NAFL-affected men, respectively. The statistics of the differential expression were calculated and the p-values were adjusted using the Benjamini & Hochberg (BH) method (12).

Next, based on the differential gene expression results, we performed a functional enrichment analysis using gene set analysis (GSA) (13). First, the genes were ordered according to their p-value and the sign of the contrast statistic. Next, the GSA was performed using the logistic regression model implemented in the mdgsa R package (14) along with their corresponding functional annotations obtained from the Gene Ontology (GO) (15) and Kyoto Encyclopaedia of Genes and Genomes (KEGG) PATHWAY (16) databases.

Because of their hierarchical structure, the gene annotations with GO terms applied in the mdgsa package were propagated so that they inherited the annotations of the ancestor terms. Excessively specific or generic annotations (blocks smaller than 10 or larger than 500 words) were subsequently filtered out. Finally, functions with a BH-adjusted p-value under 0.05 were considered significant.

For the two functional elements (GO terms and KEGG paths), we analysed the number of over-represented elements shared by the studies. These results were graphically represented as UpSet plots (17) to show the number of elements in common between the different sets. First, we compared the overrepresented elements as separate graphs to detect the common functions in each group: GO terms in premenopausal women, GO terms in men, KEGG pathways in premenopausal women, and KEGG pathways in men. We then visualised any significant GO terms in each of the individual studies in the same plot to highlight significant terms with different signs from among the studies.

5. Differential pathway activation and functional analysis

The Hipathia algorithm (18) was used to perform Pathway Activity Analysis (PAA). This method transforms gene expression values for stimulus-response signalling sub-pathways into the activation levels which ultimately trigger cellular responses. In this current work we used the hipathia R package to analyse 1,654 GO terms and 142 UniProt functions (19), associated with 146 KEGG routes (16).

The activation signal of each subpathway was computed from gene expression values using this algorithm. These values were used to detect significant differential activations in the (NASH.W – NAFL.W) – (NASH.M – NAFL.M) contrast pair (as defined above). The functional annotations included in Hipathia and the differential activation results were integrated using the GSA method to identify differences at the functional level (GO terms and UniProt functions) in relation to progression of the disease.

5. Meta-analysis

Once the gene functional analysis was applied to each individual study, we carried out a functional meta-analysis to summarise the results. Similarly, once the pathway functional analysis was applied to each individual study, the results were summarised via a homologous functional meta-analysis.

Both meta-analyses were performed following the methodology described in detail in García-García (20). Briefly, the metafor R package (21) was used to evaluate the combined effect with a random-effects model; this model more precisely detects overrepresented elements than performing the studies individually at a higher statistical power. Likewise, variability in the individual studies was incorporated into the global estimation of the measured effect so that less-variable results had a higher weight in the overall calculation of the logarithm of the odds ratio (LOR). Finally, we confirmed the suitability of the selected studies by evaluating the heterogeneity of these indicators.

For each meta-analysis and function, p-values (adjusted using the BH method), LORs, and 95% confidence intervals (CIs) were calculated. A term was considered significant when its adjusted p-value was lower than 0.05; a positive LOR in significant functions indicated a greater overrepresentation in premenopausal women than in men. In contrast, a negative LOR indicated higher representation in men than in premenopausal women. Finally, funnel plots and forest plots were used to assess the variability and to measure the contribution of each study to the meta-analysis.

A total of 8,223 elements were analysed in the gene functional meta-analysis (7,994 GO terms and 229 KEGG pathways); 1,654 GO terms and 142 UniProt functions were assessed in the functional pathway meta-analysis (20) in association with 146 KEGG pathways (16).

Fig. 1. Data-analysis workflow.

© 2020 by the authors. Submitted for possible open access publication under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).