Abstract
Introduction: Schizophrenia (SCZ), a severe neuropsychiatric disorder with high genetic susceptibility, has high rates of misdiagnosis due to the unavoidably subjective factors and heterogeneous clinical presentations. Hypoxia has been identified as an importantly risk factor that participates in the development of SCZ. Therefore, development of a hypoxia-related biomarker for SCZ diagnosis is promising. Therefore, we dedicated to develop a biomarker that could contribute to distinguishing healthy controls and SCZ patients. Methods: GSE17612, GSE21935, and GSE53987 datasets, consisting of 97 control samples and 99 SCZ samples, were involved in our study. The hypoxia score was calculated based on the single-sample gene-set enrichment analysis using the hypoxia-related differentially expressed genes to quantify the expression levels of these genes for each SCZ patient. Patients in high-score groups were defined if their hypoxia score was in the upper half of all hypoxia scores and patients in low-score groups if their hypoxia score was in the lower half. GSEA was applied to detect the functional pathway of these differently expressed genes. CIBERSORT algorithm was utilized to evaluate the tumor-infiltrating immune cells of SCZ patients. Results: In this study, we developed and validated a biomarker consisting of 12 hypoxia-related genes that could distinguish healthy controls and SCZ patients robustly. We found that the metabolism reprogramming might be activated in the patient with high hypoxia score. Finally, CIBERSORT analysis illustrated that lower composition of naive B cells and higher composition of memory B cells might be observed in low-score groups of SCZ patients. Conclusion: These findings revealed that the hypoxia-related signature was acceptable as a detector for SCZ, providing further insight into effective diagnosis and treatment strategies for SCZ.
Introduction
Schizophrenia (SCZ), as a severe and chronic psychiatric disorder, affects 1% of population worldwide [1]. The SCZ has an enormous influence on individual and their families, causing a serious socioeconomic burden. Current SCZ diagnosis is based on clinical judgment from psychiatrists, taking into account symptom presentation, structured interviews, and information from caretakers [2]. However, on account of the unavoidably subjective factors and heterogeneous clinical presentations, these diagnostic criteria usually induce the high rates of misdiagnosis [3, 4], which has fueled intensive efforts to identify a biomarker for auxiliary diagnosis of SCZ.
SCZ is a disease with high genetic susceptibility and modulated by the interplay of genetic and environmental risk factors [5]. Previous study showed that SCZ gave a life-time risk of ∼1% and showed high heritability (∼69–81%) [6, 7]. Thus, a gene biomarker derived from genome-wide expression patterns of SCZ may be promising for realizing accurate and effective diagnosis of SCZ. Lots of candidate gene biomarkers have been unveiled for SCZ [8‒10], but accurate and reliable biomarkers are still lacking.
Hypoxia has been identified as a risk factor that participates in the development of SCZ, and hypoxia-inducible factor 1 usually played a vital role [11, 12]. Candidate genes for SCZ have been revealed to subject to ischemia-hypoxia regulation or associate with vascular expression [13]. Analysis of genomic data [14] and exome sequencing [15] has shown that HIF-1A is correlated to SCZ, while pathway annotations in other genomic studies of SCZ also highlighted hypoxia genes [16]. Therefore, identifying a hypoxia-related biomarker to diagnose SCZ is beneficial to eliminate subjective factors in current diagnostic criteria and optimize SCZ diagnosis. However, no study has been reported dedicated to develop a hypoxia-related biomarker in a comprehensive, genome-wide profiling study of SCZ. Therefore, our purpose in this study is to excavate a diagnostic signature consisting of hypoxia-related genes for SCZ.
In this study, we developed a robustly hypoxia-related signature consisting of 12 genes that could distinguish healthy controls and SCZ patients in GSE17612, GSE21935, and GSE53987. Moreover, MCODE analysis indicated that vascular endothelial growth factor (VEGFA), 2-oxoglutarate 5-dioxygenase 2 (PLOD2), prolyl 4-hydroxylase subunit alpha 1 (P4HA1), and lysine demethylase 3A (KDM3A) were the remarkable nodes and highly related to SCZ as they had higher degree values in the protein-protein interaction (PPI) network of these 12 genes. Further analysis revealed that the metabolism pathway might play key roles in connecting hypoxia with the emerging of SCZ. Finally, CIBERSORT analysis illustrated that lower composition of naive B cells and higher composition of memory B cells might be observed in low-score groups of SCZ patients, suggesting that naive B cells and memory B cells might play vital roles in the pathogenesis of hypoxia-induced SCZ. Our findings provide an effective strategy to contribute to distinguishing healthy controls and SCZ patients, which not only offer a foundation for subsequent, in-depth immune-related work with great promise for the personalized treatment of SCZ but also pave the way for further investigation of mechanism for SCZ development.
Materials and Methods
Patient Data
Gene expression data and associated clinical characteristics of SCZ patients were downloaded from the publicly available GEO database (https://www.ncbi.nlm.nih.gov/geo/) (GSE17612, GSE21935, and GSE53987). Microarray gene expression profiling was GSE17612, GSE21935, GSE53987 conducted by GPL571 (Affymetrix Human Genome U133A 2.0 Array) from the GEO database. The list of these datasets is displayed in Table 1.
List of GEO datasets
Accession . | Title . | Data type of expression . |
---|---|---|
GSE17612 | Comparison of post-mortem tissue from brain BA10 region between schizophrenic and control patients | Microarray data |
GSE21935 | Comparison of post-mortem tissue from Brodmann brain BA22 region between schizophrenic and control patients | Microarray data |
GSE53987 | Microarray profiling of PFC, HPC, and STR from subjects with SCZ, bipolar, MDD, or control | Microarray data |
Accession . | Title . | Data type of expression . |
---|---|---|
GSE17612 | Comparison of post-mortem tissue from brain BA10 region between schizophrenic and control patients | Microarray data |
GSE21935 | Comparison of post-mortem tissue from Brodmann brain BA22 region between schizophrenic and control patients | Microarray data |
GSE53987 | Microarray profiling of PFC, HPC, and STR from subjects with SCZ, bipolar, MDD, or control | Microarray data |
Differentially Expressed Genes Related to Hypoxia Screening in GSE17612
First, we performed differential gene analysis of 21,655 mRNAs using the package “limma” for comparisons of SCZ samples and controls in GSE17612. |LogFC| > 2 and false discovery rate (FDR) < 0.05 were set as the cut-off values. Next, Homo sapiens hypoxia-related gene set was defined based on Chemical and Genetic Perturbations (CGP) (MANALO_HYPOXIA_UP), which included 206 genes. These gene sets were presented in online supplementary Table S2 (see online suppl. material at www.karger.com/doi/10.1159/000529902). Then, we selected common genes which were overlapped with differentially expressed genes (DEGs) and list of hypoxia-related genes for further analysis.
Construction of Hypoxia-Related Signature Based on Hypoxia Score
The hypoxia score was calculated based on the single-sample gene-set enrichment analysis (ssGSEA) [17] using the 12 hypoxia-related DEGs (online suppl. Table S3) to quantify the expression levels of these genes for each sample. We estimated the hypoxia score between SCZ and control samples in GSE17612. Additionally, SCZ patients as high-score group were defined if their hypoxia score was in the upper half of all hypoxia scores and others as low-score group if their hypoxia score was in the lower half.
Protein-Protein Interaction Network and Functional and Pathway Enrichment Analysis
A PPI network of the hypoxia-related DGEs was constructed by the Search Tool for the Retrieval of Interacting Genes database (https://string-db.org) and visualized by Cytoscape software. For a more in-depth investigation of the functions of the hypoxia-related signature, the Database for Annotation, Visualization, and Integrated Discovery database (https://david.ncifcrf.gov/) was implemented to enrich biological themes on gene ontology (GO) terms.
To explore the potential biological functions of the hypoxia-related gene set in SCZ development, we investigated the biological pathway enriched in high-score groups using gene-set enrichment analysis (GSEA) algorithm. SCZ patients were divided into “high” and “low” hypoxia groups, which were determined by the cut-off value of median hypoxia score based on the expression of 12 hypoxia-related DEGs. Then, all of SCZ patients in “high” and “low” hypoxia groups were put into GSEA analysis based on the curated gene sets c2.cp.kegg.v7.0.symbols using GSEA version 4.1.0 software package. Enrichment score reflects the enrichment degree of a pathway-related gene set. The larger the enrichment score, the more relevant the reaction is to this pathway. p value < 0.05 was considered to be statistical significance.
Evaluation of SCZ-Infiltrating Immune Cells
Normalized gene expression data were used to evaluate the relative proportions of 22 types of infiltrating immune cells using the CIBERSORT algorithm [18]. Gene expression datasets were prepared using standard annotation files and data uploaded to the CIBERSORT Web portal (https://cibersort.stanford.edu), with the algorithm run using the default signature matrix at 1,000 permutations, which has been applied in many studies [19‒22]. Here, we first divided SCZ patients into “high” and “low” hypoxia groups, which were determined by the cut-off value of median hypoxia score based on the expression of 12 hypoxia-related DEGs. And then, we applied the original CIBERSORT gene signature file LM22 which defines 22 immune cell subtypes and analyzed the immune cell composition in these SCZ samples from GSE17612 dataset.
Statistical Analysis
The expression profiles of the mRNAs from GEO are shown as the raw data; then each mRNA was normalized by log2 transformation for further analysis. The receiver operating characteristic (ROC) curves were plotted based on the risk scores and the case-control classification of each patient to compare the predictive accuracy of the gene signature. p values from the log-rank tests were calculated, and a p value < 0.05 was considered statistically significant. Statistical analysis was performed by using GraphPad Prism version 7.0 or SPSS version 19.0 software package. A two-tailed p < 0.05 or FDR <0.05 was considered statistically significant.
Results
Identification and Functional Analysis of Differentially Expressed Genes in SCZ and Healthy Control Samples
We utilized GSE17612 from GEO dataset to identify DEGs in 28 SCZ patients and 23 healthy controls. After screening using the criteria of |logFC| > 2 and FDR <0.05, 885 genes with significantly different expression levels between the two groups (466 up-regulated genes and 419 down-regulated genes) were unveiled (online suppl. Table S1; Fig. 1a). Subsequently, functional enrichment analysis revealed that these genes were most enriched in several GO terms related to energy metabolism. Specifically, GO analysis (biological processes, cellular components, and molecular functions) showed that the DEGs were mainly enriched in “oxidation-reduction process,” “mitochondrion,” and “protein homodimerization activity,” respectively (Fig. 1b, c, d). For the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, the DEGs were mainly concentrated in metabolic pathway (Fig. 1e).
Functional enrichment analysis of differentially expressed genes (DEGs) in GSE17612. a Volcano plot indicated up-regulated and down-regulated mRNAs in SCZ patients. Red represents the selected up-regulated genes and green represents the selected down-regulated genes. Gene ontology (GO) analysis for DEGs: biological processes (BPs), molecular function (MFs), and cellular components (CCs) (b–d). e Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enrichment of DEGs.
Functional enrichment analysis of differentially expressed genes (DEGs) in GSE17612. a Volcano plot indicated up-regulated and down-regulated mRNAs in SCZ patients. Red represents the selected up-regulated genes and green represents the selected down-regulated genes. Gene ontology (GO) analysis for DEGs: biological processes (BPs), molecular function (MFs), and cellular components (CCs) (b–d). e Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enrichment of DEGs.
Identification of DEGs Related to Hypoxia
As we learned in the introduction, hypoxia might be an important cause of SCZ [11, 12]. Therefore, it is promising to excavate a hypoxia-related signature for distinguishing SCZ from controls. Homo sapiens hypoxia-related gene set was defined based on Chemical and Genetic Perturbations (CGP: MANALO_HYPOXIA_UP), which included 206 genes. Out of 885 DEGs, we extracted 12 genes that overlapped with the genes in the list of hypoxia genes. These 12 common genes included 3 up-regulated and 9 down-regulated genes (Fig. 2a, b).
Identification of hypoxia-related DEGs in GSE17612. a Venn diagram displayed 12 hypoxia-related genes overlapped with DEGs. b Heat map showed the different expression levels of the 12 hypoxia-related genes between SCZ samples and controls. c Gene ontology (GO) analysis for hypoxia-related DEGs: biological processes (BPs), molecular function (MFs), and cellular components (CCs). d PPI network of hypoxia-related DEGs constructed by the STRING database. e Identification of a sub-network using MCODE analysis in Cytoscape software.
Identification of hypoxia-related DEGs in GSE17612. a Venn diagram displayed 12 hypoxia-related genes overlapped with DEGs. b Heat map showed the different expression levels of the 12 hypoxia-related genes between SCZ samples and controls. c Gene ontology (GO) analysis for hypoxia-related DEGs: biological processes (BPs), molecular function (MFs), and cellular components (CCs). d PPI network of hypoxia-related DEGs constructed by the STRING database. e Identification of a sub-network using MCODE analysis in Cytoscape software.
Subsequently, the 12 genes in the biological process group were mainly enriched in “oxidation-reduction process” and “response to hypoxia”; in the CC group, the genes were mainly concentrated in “endoplasmic reticulum”; in the molecular function group, the genes were mainly enriched in “iron-ion binding” and “l-ascorbic acid binding” (Fig. 2c). These results confirmed that our genes were closely associated with hypoxia. To better understand the interplay among these hypoxia-related DEGs, we constructed a PPI network (Fig. 2d) and screened out only one module with score value of four in MCODE analysis, including VEGFA, PLOD2, P4HA1, and KDM3A. These genes were the remarkable nodes and could be seen as hub genes as they had higher degree values in this network (Fig. 2e).
The Hypoxia-Related Signature Has Potential to Diagnose SCZ
We performed ssGSEA to calculate the hypoxia score based on the 12 genes expression data and compared the hypoxia score between SCZ and control samples using GSE17612 dataset (Fig. 3a). The results showed that SCZ samples had significantly different hypoxia scores compared to control samples in GSE17612 dataset (t = 2.836, p = 0.0066, Fig. 3b). To validate these results, we also involved in 2 independent datasets, and the results were consistent with those in GSE17612 dataset (GSE21935: t = 2.771, p = 0.0084; GSE53987: t = 3.607, p = 0.0005, Fig. 3a, b), which strengthened the generalizability of our findings.
Construction and validation of a hypoxia-related signature for SCZ diagnosis in the GSE17612, GSE21935, and GSE53987 datasets. a Significant difference of hypoxia score between SCZ and control samples in GSE17612, GSE21935, and GSE53987 datasets. * p < 0.05; ** p < 0.01; *** p < 0.001. b Distribution of patient hypoxia scores for SCZ patients and controls. c ROC curve analysis reveals the predictive value of the hypoxia-related signature for SCZ patients in the GSE17612, GSE21935, and GSE53987 datasets.
Construction and validation of a hypoxia-related signature for SCZ diagnosis in the GSE17612, GSE21935, and GSE53987 datasets. a Significant difference of hypoxia score between SCZ and control samples in GSE17612, GSE21935, and GSE53987 datasets. * p < 0.05; ** p < 0.01; *** p < 0.001. b Distribution of patient hypoxia scores for SCZ patients and controls. c ROC curve analysis reveals the predictive value of the hypoxia-related signature for SCZ patients in the GSE17612, GSE21935, and GSE53987 datasets.
Given the significant transcriptional differences between SCZ patients and healthy controls, we hypothesized that the hypoxia score based on the 12 DEGs could be generated for SCZ detection. Therefore, we tested this hypothesis using ROC curves in GSE17612, GSE21935, and GSE53987 datasets. Interestingly, the results showed that the hypoxia-related signature based on the 12 DEGs achieved AUC of 0.763, 0.713, and 0.636 for GSE17612, GSE21935, and GSE53987, respectively, in diagnosing SCZ, suggesting moderate potential of the hypoxia score as a detector of SCZ (Fig. 3c). Above all, these data indicated that the hypoxia score based on the 12 DEGs had potential for recognizing the SCZ patients.
The Hypoxia-Related Signature Is Related to Metabolism Pathways in SCZ
To investigate the underlying mechanism of the hypoxia score in diagnosing SCZ patients, we compared the gene expression profiles between high- and low-score groups and found that the metabolism-related pathways were up-regulated in high-score groups, including arginine and proline metabolism, lysine degradation, and glycolysis gluconeogenesis (Fig. 4a, b). Altogether, these data indicate that some of these metabolism-related pathways may play key roles in connecting hypoxia with incidence of SCZ. SCZ patients with high hypoxia score are prone to have high composition of naive B cells and low composition of memory B cells.
Genes expressed in high-score groups are mostly enriched in metabolism-related BPs. a List of the top 20 KEGG pathways using GSEA analysis. b KEGG pathways of the hypoxia-related signature in SCZ patients using GSEA analysis.
Genes expressed in high-score groups are mostly enriched in metabolism-related BPs. a List of the top 20 KEGG pathways using GSEA analysis. b KEGG pathways of the hypoxia-related signature in SCZ patients using GSEA analysis.
A mass of researches have proved that hypoxia usually affected immune functions [23, 24]. Given this, one emerging question is whether the hypoxia signature is associated with immune cell composition of SCZ patients. Therefore, we performed further analysis using CIBERSORT algorithm to illuminate the relationships between the hypoxia score and immune cell composition using GSE17612. Significantly, each SCZ patient had various proportions of immune cells in high- and low-score groups (Fig. 5a, b). According to the result of Figure 5c, the proportions of infiltrating immune cells subpopulations were weakly to moderately corelated in SCZ patients. Surprisingly, as shown in Figure 5d, naive B cells showed significantly positive correlation with the hypoxia scores, while memory B cells showed significantly negative. These results demonstrated that lower composition of naive B cells and higher composition of memory B cells might be observed in low-score groups of SCZ patients.
a Summary of immune composition of the immune cells’ subpopulations in SCZ samples in GSE17612. b Heat map of the immune cell proportions in GSE17612. c Correlation matrix of the immune cell proportions in GSE17612. d Different fractions of the immune cells in low- and high-score groups in GSE17612.
a Summary of immune composition of the immune cells’ subpopulations in SCZ samples in GSE17612. b Heat map of the immune cell proportions in GSE17612. c Correlation matrix of the immune cell proportions in GSE17612. d Different fractions of the immune cells in low- and high-score groups in GSE17612.
Discussion
Currently, the classic diagnostic method is clinical interview and WHO-DAS score which are widely used clinically [2, 25]. But these two methods were usually affected by subjective factors, leading to decrease the authenticity and reliability of diagnosis. Thereby, more and more studies paid attention to seek novel biomarkers for screening SCZ. Despite a few studies had made progress, there were still some limitations which reduced the reliability of conclusions. SCZ is a severe psychiatric disorder with a genetic susceptibility [26]. Moreover, hypoxia has been recognized as a risk factor for increase of SCZ incidence [13, 15, 27]. Given these data, a hypoxia-related biomarker is promising for SCZ diagnosis. Therefore, in this study, we dedicated to identify a gene biomarker related to hypoxia for SCZ diagnosis based on gene expression profiles from GEO database.
We screened out 885 differential genes and analyzed these DEGs through a variety of bioinformatics analysis methods. GO terminology indicate that these 885 DEGs were enriched in biological processes, such as “oxidation-reduction process,” “mitochondrion,” and “protein homodimerization activity.” Previous study implicated redox dysregulation as a pathological mechanism driving the emergence of SCZ [28]. Researchers indicated that increased oxidative damage and decreased capacity of intracellular redox modulatory systems were consistent findings in persons with SCZ as well as in persons at clinical high risk who subsequently developed frank psychosis [29]. Multiple pieces of evidence suggested that brain energy metabolism, mitochondrial functions, and redox balance were impaired to various degrees in SCZ [30]. KEGG pathway enrichment analysis shows that these DEGs are mainly related to “metabolic pathway.” Our findings suggest that the “oxidation-reduction process,” “mitochondrion,” “protein homodimerization activity,” and “metabolic pathway” may play a vital role in the pathogenesis of SCZ.
Next, in view of the vital role in driving the emergence of SCZ, we finally excavated 12 DEGs which is related to hypoxia pathway. Thus, a hypoxia-related signature comprising12 DEGs was established. Subsequently, functional enrichment analysis revealed that these genes were enriched in several GO terms related to “oxidation-reduction process,” “response to hypoxia,” “endoplasmic reticulum,” “iron-ion binding,” and “l-ascorbic acid binding,” suggesting that our genes were indeed associated with hypoxia pathway. A PPI network was further implemented to study the 12 genes’ interaction relationship and screened out the modules with the highest scores in MCODE analysis, including VEGFA, PLOD2, P4HA1, and KDM3A. These genes were the remarkable nodes and highly related to SCZ as they had higher degree values in this network. Consistently, previous study unveiled that VEGFA showed a significantly higher expression variance among SCZ samples than control samples [31]. Paulo Lizano et al. [32] reported that VEGFA dysfunction might contribute to a number of pathological processes that characterize SCZ. Procollagen lysine, PLOD2 have been confirmed to have an anti-metastasis and anti-proliferation effect [33]. P4HA1 has been proved to associate with poor prognosis and promote cell proliferation and metastasis of adenocarcinoma [34, 35]. Previous study proved that KDM3A could protect the myocardium from ischemia/reperfusion injury [36] and promote osteosarcoma development [37]. However, in addition to VEGFA, other three genes (PLOD2, P4HA1, and KDM3A) have not been reported in SCZ. To the best of our knowledge, this is the first study to unveil the potential power of PLOD2, P4HA1, and KDM3A as HIF-1α-related genes for the SCZ development. This new finding points a new direction for SCZ research in the future.
Subsequently, ssGSEA was employed to calculate the hypoxia score based on the 12 genes expression data, and we compared the hypoxia score between SCZ and control samples using GSE17612, GSE21935, and GSE53987 datasets. Surprisingly, the hypoxia-related signature displayed significant difference between SCZ and control samples, which indicated potential value for SCZ detection. Then, ROC curves confirmed that the hypoxia-related signature based on the 12 DEGs had moderate potential as a detector for SCZ.
To investigate the underlying mechanism of the hypoxia-related signature in diagnosing SCZ, we further carried out KEGG pathway analysis. The results showed that the metabolism-related pathways were up-regulated in high-score groups, including arginine and proline metabolism, lysine degradation, and glycolysis gluconeogenesis. Altogether, these data indicate that some of these metabolism-related pathways may play key roles in connecting hypoxia with the emerging of SCZ.
Hypoxia is now recognized to contribute fundamentally to immune responses. A large number of evidence have implicated a prominent role for hypoxia in innate immunity of healthy tissue (physiologic hypoxia) and during active inflammation (inflammatory hypoxia) [38, 39]. Given this, we reasoned that if the hypoxia-related signature related to immune responses, we set out to test this hypothesis by performing CIBERSORT algorithm to illuminate the relationships between the hypoxia score and immune cell composition using GSE17612. Significantly, naive B cells showed significantly positive correlation with the hypoxia scores, while memory B cells showed significantly negative. These results demonstrated that lower composition of naive B cells and higher proportion of memory B cells might be observed in low-score groups of SCZ patients, suggesting that naive B cells and memory B cells might play vital roles in the pathogenesis of hypoxia-induced SCZ.
Despite the significant results obtained in the current study, there were inevitably several shortcomings of our study that should be acknowledged. First, gene expression data of SCZ patients were downloaded from publicly available datasets. However, the information of publicly available datasets is limited, so that the demographic variables could be potential confounders, which might induce bias results. Second, we did not consider the heterogeneity of the immune microenvironment related to the location of immune infiltration. Third, there were no experimental data regarding the identified signature. Therefore, further research is needed to elucidate the inherent correlation between the 12 hypoxia-related genes and the development of SCZ.
Conclusion
In conclusion, we developed and validated a hypoxia-related signature consisting of 12 DEGs that could contribute to distinguishing between healthy controls and SCZ patients. Metabolism reprogramming might be activated in the patient with high hypoxia score. Additionally, patients in high hypoxia-score groups had more composition of naive B cells and lower proportion of memory B cells.
Statement of Ethics
The original TCGA data that support the findings of our study are available in the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) data from NCBI GEO (accession numbers: GSE17612, GSE21935, and GSE53987). An ethics statement is not applicable because this study is based exclusively on published literature.
Conflict of Interest Statement
The authors declare that they have no conflict of interest.
Funding Sources
There is no funding to report for this manuscript.
Author Contributions
Conceptualization: Zhitao Li; collection and assembly of data: Xinyu Sun; data analysis and interpretation: Jia He; original draft preparation: Dongyan Kong; writing: Jinyi Wang; and visualization and supervision: Lili Wang.
Data Availability Statement
The original TCGA data that support the findings of our study are available in the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) data from NCBI GEO (accession numbers: GSE17612, GSE21935, and GSE53987). Further inquiries can be directed to the corresponding author.