Microarray-based gene expression profile studies were identified in the Pubmed database (http://www.pubmed.com), Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/gds/) and ArrayExpress dataset of the European Molecular Biology Laboratory-European Bioinformatics Institute (http://www.ebi.ac.uk/arrayexpress/) for MDS and AML. The following key words were used for MDS and AML respectively: myelodysplastic syndrome or dysmyelopoietic syndrome or hematopoetic myelodysplasia or MDS, and microarray or gene expression profile or gene expression profiling and acute myeloid leukemia or acute myeloblastic leukemia or acute myelocytic leukemia or acute nonlymphoblastic leukemia or acute nonlymphocytic leukemia) or acute myelogenous leukemia or AML or ANLL and microarray or gene expression profile or gene expression profiling.

Eligible studies and datasets should follow these inclusive criteria: i) patient and healthy control studies of human; ii) analysis of gene expression profiling; iii) comparable experimental conditions and untreated; and iv) available complete raw and processed microarray data. Studies were excluded if they were: i) letters, abstracts, meta-analysis, review articles and case repor; ii) cell lines were used in experimental design; iii) RT-PCR only for profiling studies; and iv) studies without healthy control. All the datasets and references, which are conformed to the criteria mentioned above, were carefully screened. The latest search was performed on May 30, 2017.

Full text and supplementary materials of selected articles were screened and the key items were extracts as listed: GEO series accession number, type of disease, number of patients and healthy donors, specimen sources and platform of microarray (Table I). The series matrix files were downloaded from GEO datasets for all studies, with the exception of GSE983, whose CEL files were obtained and processed further by R platform to generate preliminary series matrix file. Common Entrez IDs were used to substitute all the gene probes in accordance with the corresponding microarray platforms. Before integrative meta-analysis, individual dataset was normalized by log2 transformation and R-mediated mean, and quantile normalization. Expression data of patients with MDS or AML and healthy donors were defined as class 2, and 1, respectively, according to the guidelines of NetworkAnalyst (10).

NetworkAnalyst is capable of massive datasets integration on the premise of batch effect adjust option. The processed and normalized datasets were uploaded and subjected to the well-established ComBat procedures to remove study-specific batch effects, which uses the Emperical Bayes method to adjust the extreme expression ratios, alleviate gene variances across all other genes, and possibly removing their inference without compromising the biological covariates (12). The sample clustering patterns with and without batch effect adjustment were visualized and compared by principal component analysis (PCA) to assess the efficiency of batch effect removal.

We conducted the meta-analysis using NetworkAnalyst, a web interface for integrative statistical and visualizing tool. In the option of ‘multiple gene expression data’ for the web interface, All datasets were uploaded to the ‘multiple gene expression data’ input area and analyzed in a streamlined manner, including data processing for Entrez ID, annotation check, view by both Boxplot and PCA plot, confirmation of normalization, individual DEGs analysis and data summary. For the DEGs discrimination, the cut-off of P-value was adjusted to 0.05, using the false discovery rate (FDR) based on Benjamini-Hochberg procedure and moderated t-test based on the Limma algorithm. Furthermore, all datasets were subjected to integrity check to ensure that the merged data could be carried out by ES combination, which allowed the generation of more biologically consistent meta-based DEGs by incorporation of both the magnitude and direction of gene fold change. Between the two popular methods: Fixed (FEM) and random effects models, REM was used in the current study. REM model assumes that each study contains a random ES that incorporates unknown cross-study heterogeneities, as demonstrated by Cochran's Q tests. ‘Define custom signature’ tool from NetworkAnalyst was used to produce the heatmap visualization for both top 25 up- and downregulated genes.

PPI Network-based analysis was performed using NetworkAnalyst (10). In PPI networks, nodes stand for proteins, while edges represent known interactions between the linked proteins. Topology and subnetwork analyses were performed in this study, in purpose to demonstrate its overall structural properties and highlight part of the network that had shown significant changes. In short, the subnetwork analysis was executed in three steps: i) obtaining the list of DEGs by meta-analysis; ii) proceeding the list to the IMEx interactom-based PPI network analysis; and iii) choosing zero-order and minimum network to avoid ‘hairball effect’ and overall presentation of the structure. In addition, there are two kinds of complementary measurements in NetworkAnalyst to reveal the most important nodes, also called hub genes: Degree and betweenness centrality. Degree centrality is the number of connections that a node has to other nodes, whereas betweenness centrality corresponds to the number of shortest paths passing through the node. From the parent network, the most significant modules of hub genes were extracted for both up- and downregulated DEGs using the ‘module explorer tool’, based on the random walks-dependent Walktrap algorithm.

To demonstrate the implication of shared DEGs in MDS and AML, we conducted an enriched pathway group analysis by ClueGO of Cytoscape, a software with significantly expanded annotated gene sets of Kyoto Encyclopedia of Genes and Genomes (KEGG) and Reactome pathway databases. The processes were submitted to the followed settings and parameters: enrichment (right-sided) hyper-geometric distribution tests, with a P-value ≤0.05; Benjamini and Hochberg adjustment for the terms and the groups with Kappa-statistics score threshold set to 0.4; leading term groups were selected based on the highest significance. The gene ontology (GO) terms enrichment analysis was conducted by topGO of Bioconductor version 2.14, R 3.1.1, and summarized by REVIGO. The enriched GO terms or pathways were considered to be significant with adjusted P-value <0.05.

The web-server NetworkAnalyst analyzed the data and performed meta-analysis by REM-based combined ES statistics. DEGs were selected by FDR guided Benjamini-Hochberg procedure with adjusted P-value <0.05. Hypergeometric test (right-sided) and Benjamini & Hochberg FDR correction were used to identify significantly enriched biological pathways in DEGs by ClueGO of Cytoscape. Significantly enriched biological processes were monitored using top GO and REVIGO analysis and P<0.05 was considered to indicate a statistically significant difference.

2012 and 476 literatures were retrieved from Pubmed for AML and MDS, respectively, according to our search strategy. And after a series of inclusion and exclusion screening, a total of 8 and 10 microarray datasets were selected for MDS (GSE19429, GSE2779, GSE30195, GSE41130, GSE4619, GSE51757, GSE58831 and GSE81173) (8,13–17) and AML (GSE24395, GSE30029, GSE33223, GSE37307, GSE68172, GSE70284, GSE79605, GSE84881, GSE9476 and GSE983) (18–25) respectively, as indicated in Materials and methods (Fig. 1A). All the 18 datasets were carefully screened to meet the inclusion and exclusion criteria for meta-analysis. Collectively, a total of 106/281 and 75/465 (healthy control/patient) samples were included for MDS and AML, respectively. There are 8 different microarray platforms applied in the selected datasets, which include Affymetrix Human Genome U133 Plus 2.0 Array, Affymetrix Human Genome U133A Array, Affymetrix Human Full Length HuGeneFL Array, Affymetrix Human Gene Expression Array, Illumina HumanHT-12 v.3.0 expression beadchip, Sentrix Human-6 v2 Expression BeadChip, Agilent-014850 Whole Human Genome Microarray, and Agilent-028004 SurePrint G3 Human GE 8×60 K Microarray. In addition, various tissues, such as total BM cells, FACS-sorted CD34⁺ hematopoietic progenitors, and peripheral blood cells, were used for the microarray analysis. In the process of data integration, patient samples of MDS and AML were not distinguished, for the purpose to uncover the common gene signature shared between MDS and AML. Table I presents the detailed information for each dataset, including GEO accession number, type of disease, sample composition, sample source, and the corresponding references.

The workflow for meta-analysis used in the present study was illustrated in Fig. 1B. To identify the shared transcriptional signatures between MDS and AML, the total of 18 datasets were simultaneously analyzed by NetworkAnalyst. When the 18 microarray datasets were analyzed individually, a total of 7579 DEGs were revealed. For meta-analysis, we uploaded the 18 datasets in succession, and each dataset was processed by Entrez ID matching, sample (control/patient) annotation and individual DEGs identification. To remove batch effects among different datasets, we performed ‘ComBat’-based batch effect adjustment (12), and the sample clustering patterns with or without batch effect adjustment were visualized by PCA plots (Fig. 2A and B). All the 18 gene expressional microarray data were then integrated and merged. Meta-analysis was conducted following the Cochran's Q test, REM and ES statistical methods, which facilitated to reveal the DEGs between healthy donors and patients across different microarray datasets by permitting variable true ES and integrating unknown cross-study heterogeneities. Finally, we found a total of 330 DEGs, including 191 up- and 139 downregulated genes across the 18 datasets with significance threshold of adjusted P-value <0.05. By comparing the list of individually identified DEGs with meta-analysis based DEGs, 211 DEGs were revealed by both analyses. Of note, 119 DEGs were uniquely discovered by meta-analysis, which were referred as gained genes, and 7378 DEGs were only found in individual analysis, which were defined as lost genes (Fig. 2C). Expression profile of the top 25 up- and downregulated genes among the 330 identified DEGs was visualized by heatmap (Fig. 2D). Due to the large number of healthy donor and patient samples, the heatmap was divided into 3 parts by GEO datasets. The gene expressional profiles of the top 25 up- and downregulated genes could be easily visualized to be consistent across different datasets, except those with small sample numbers, including GSE33223, GSE79605, GSE41130 and GSE84881. The Parathyroid hormone 2 receptor (PTH2R), Tec protein tyrosine kinase (TEC), and Glutathione peroxidase 1 (GPX1) were among the most significantly upregulated genes, and PTH2R had the highest combined ES of 1.0994, which was upregulated consistently with ES ranging from 0.47068 to 2.8508 across all the 18 datasets (data not shown), while Matrix metalloproteinase 12 (MME), Recombination activating 1 (RAG1), and CD79b molecule (CD79B) were among the most significantly downregulated genes, and MME had the relatively highest combined ES of −1.29, which was consistent in 15 out of 18 datasets with ER ranging from −3.7844 to −0.33961 (data not shown). Additionally, top 10 up- and downregulated DEGs were listed in Table II.

Many studies have shown that ‘hub’ proteins are more likely to be encoded by pleiotropic genes or genes that are related to certain diseases (26). Hub nodes are potentially key molecules in signaling, as they are highly interconnected with dysregulated genes and they often receive and integrate multiple signals and pass them onto downstream nodes (27). In order to identify the key hub genes among the common DEGs shared by MDS and AML, we intentionally performed a network based meta-analysis. NetworkAnalyst provided us with a protein-protein interaction (PPI) network analysis tool by integrating the IMEx (International Molecular Exchange Consortium) interactome with the original seed of 330 DEGs. The NetworkAnalyst-based PPI network analysis generated first-order, minimum-order and zero-order networks with (4507 nodes, 9421 edges), (873 nodes, 3510 edges), and (123 nodes, 157 edges), respectively. The minimum-order network shows the overall structure (Fig. 3A), while the zero-order network helps to visualize the detailed interactions between different seeds (Fig. 3B). HSP90AA1 (combined ES: 0.36197, adjusted P-value: 0.0026418) and CUL1 (combined ES: −0.54806, adjusted P-value: 0.002118) were found to be the top ranked hub genes by degree and betweenness centrality analysis among the up- and downregulated DEGs. We also listed the top 10 ranked hub genes with detailed information in Table III. Additionally, we used the ‘module explorer’ tool to highlight two most significant modules as two sub-networks composed of HSP90AA1 (19 nodes and 23 edges) and CUL1 (20 nodes and 25 edges), respectively (Fig. 3C and D).

In order to thoroughly understand the pathways involved in MDS and AML development, the enriched biological pathways for up- and downregulated DEGs were functionally grouped with the threshold of P-value <0.05 by the ClueGO plugin of Cytoscape v.3.4.0, and only KEGG and Reactome pathway databases were selected. For upregulated DEGs, we found 88 enriched biological pathways, which can be divided into 16 groups, mainly involved in ‘Signaling by fibroblast growth factor receptor (FGFR) in disease’, ‘EPH-ephrin mediated repulsion of cells’, ‘Histidine, lysine, phenylalanine, tyrosine, proline and tryptophan catabolism’, ‘Positive epigenetic regulation of rRNA expression’, and ‘Neurophilin interactions with VEGF and VEGFR’ (Fig. 4A). For downregulated DEGs, we found 37 enriched biological pathways, which can be divided into 17 groups, mainly involved in ‘Interleukine-6 signaling’, ‘Integration of provirus’, ‘Graft-versus-host disease’, ‘B cell receptor signaling pathway’, and ‘Effects of PIP2 hydrolysis’ (Fig. 4B).

Gene set enrichment analyses were conducted for the analysis of overrepresented gene ontology (GO) terms, which were enriched by topGO package of R for the 191 up- and 139 downregulated DEGs with the threshold of P-value <0.05. REVIGO is a web server that summarizes long, unintelligible lists of GO terms by finding a representative subset of the terms using a simple clustering algorithm that relies on semantic similarity measures (28). Using REVIGO, the top 200 enriched GO terms for upregulated DEGs were shown in 20 subsets, including ‘positive regulation of tyrosine phosphorylation of Stat4 protein’ (comprised by 58 terms), ‘adrenal cortex formation’ (comprised by 33 terms), ‘establishment or maintenance of transmembrane electrochemical gradient’ (comprised by 20 terms), ‘heart contraction’ (comprised by 10 terms), etc, which were visualized by treemap generated by R (Fig. 5A). Similarly, 13 subsets were shown for downregulated DEGs, including ‘B cell receptor signaling pathway’ (comprised by 75 terms), ‘embryo implantation’ (comprised by 24 terms), ‘fructose metabolism’ (comprised by 23 terms), ‘response to lipopolysaccharide’ (comprised by 18 terms), etc, which were also visualized by treemap (Fig. 5B).