Mining database for the therapeutic targets and prognostic biomarkers among STAT family in glioblastoma

2.1 Oncomine

Oncomine (www.oncomine.org/) is a web application offering translational bioinformatics services [11]. It contains over 700 independent datasets, enabling researchers to interrogate gene expression profiles and biological interactions. Using Oncomine, mRNA levels of STATs in GBM were determined at 0.05 and a fold-change = 2.

2.2 GEPIA

Gene Expression Profiling Interactive Analysis (GEPIA – www.gepia.cancer-pku.cn/index.html) is web application developed by Zefang et al. at Peking University. Using GEPIA, customizable functions including tumor/normal differential expression analysis and correlation analysis of STAT family in GBM were explored [12]. Pearson correlation coefficient was utilized to performed correlation analysis using The Cancer Genome Atlas (TCGA), GBM dataset ( N = 163) with a p -value of 0.05. Survival curves were drawn using Kaplan-Meier method.

2.3 PrognoScan

PrognoScan (http://dna00.bio.kyutech.ac.jp/PrognoScan/index.html) is a public database for performing meta-analysis of prognostic value of genes [13]. By searching large, publicly available cancer microarray datasets, it enables evaluation of relationships between gene expression and prognosis. PrognoScan was used to explore the significance of STAT family in prognosis of GBM patients.

2.4 TIMER

Tumor IMmune Estimation Response (TIMER, www.cistrome.shinyapps.io/timer/) is web application developed by Li et al. for analyzing immune infiltrates [14]. The web server was used to determine immune cell infiltrates and immune gene biomarker correlation of STAT family in GBM. Table 3 shows immune gene biomarkers and their corresponding immune cells, which have been identified in previous studies [15, 16, 17]. The analysis was conducted using Spearman correlation coefficient based on TCGA GBM dataset ( N = 163).

2.5 GSCALite

Gene Set Cancer Analysis (GSCALite, www.bioinfo. life.hust.edu.cn/web/GSCALite/) is a web-based platform developed by Liu et al. for analyzing genes related to cancer [18]. Customizable functions including SNV, cancer pathway activity and drug analysis of STAT family in GBM were explored using GSCALite. By analyzing cancer pathway activity, STAT family members were grouped into two (High and Low) by median expression. The difference in pathway activity score (PAS) between the groups was defined using Student’s t -test and P-value adjusted by FDR. FDR 0.05 was considered significant. When PAS (Gene A group High) $>$ PAS (Gene A group Low), we considered gene A to have an activating effect on a pathway, otherwise it would have an inhibiting effect. Spearman correlation test was used to analyze correlation between STAT family members and 265 small molecules from Genomics of Drug Sensitivity in Cancer (GDSC). All the analyses were conducted using TCGA GBM dataset ( N = 163).

OPEN IN VIEWER

LinkedOmics (www.linkedomics.org) is a publicly available portal containing multi-omics data and clinical data of 11158 patients from the TCGA project, which could performed translational bioinformatics services [19]. Various customizable functions including correlated gene analysis, Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, kinase and transcription factor-target analysis were explored using LinkedOmics and the TCGA GBM sample ( N = 163). Pearson Correlation and Gene Set Enrichment Analysis (GSEA) were conducted at p -value 0.05.

2.7 GeneMANIA

GeneMANIA (www.genemania.org) is a prediction server for gene prioritization and predicting gene function [20]. Customizable functions including kinase LCK and transcription factor ELF1 network were explored with GeneMANIA by constructing a protein-protein interaction (PPI) network.

3.1 Expression and prognosis analysis of STAT family in GBM

To assess the potential involvement of STAT family in GBM, their differential expression at mRNA and protein level was analyzed using multiple datasets hosted on Oncomine and GEPIA. Analysis of STAT family mRNA expression in 20 cancer types (Oncomine) showed a significant upregulation and downregulation in several cancers including GBM (Fig. 1). Particularly, STAT1 was significantly upregulated in tumor tissues, with a fold change of 3.136 and a p -value of 1.58E-10 (Fig. 1, Table 1) [21]. This finding is consistent with a previous study that reported upregulation of STAT1 in tumor tissues with a fold change of 2.701, and a p -value of 1.69E-5 [22]. Two datasets showed an increase of STAT3 in GBM with a fold change of 2.076 and 2.270 [23, 24]. Data obtained from Lixin et al. showed STAT5A and STAT5B to be increased in GBM (STAT5A: fold change = 2.255; STAT5B fold change = 2.953) [24]. Expression of STAT6 was significantly upregulated GBM, with a fold change 2.103 ( p = 0.004). Differential expression analysis using GEPIA showed high levels of STAT1, STAT3, and STAT5A, and a low level of STAT4 in GBM tissues compared with brain tissues (Fig. 2A).

OPEN IN VIEWER

OPEN IN VIEWER

The results of prognosis analysis of STAT family in GBM using PrognoScan are shown in Table 2. Notably, GBM patients with high expression of STAT1 (HR [95% CI]: 1.55 [1.11–2.16], p = 0.010313), and STAT2 (HR [95% CI]: 1.55 [1.11–2.16], p = 0.018909) had the worse overall survival. Two datasets showed GBM patients with high expression of STAT4 (HR [95% CI]: 0.25 [0.09–0.70], p = 0.007913), and STAT5B (HR [95% CI]: 0.74 [0.46–1.19], p = 0.012642) to have a better overall survival. Data obtained from TGGA were used to verify these results, indicating that GBM patients with high expression of STAT1 ( p = 0.00035), STAT2 ( p = 0.00032), STAT3 ( p = 0.016), STAT5A ( 0.0001) and STAT6 ( p = 0.0086) had the worse overall survival (Fig. 2B).

OPEN IN VIEWER

Missense mutation was the most common variant classification and C $>$ T was the most common SNV class (Fig. 3A). Among the members of the STAT family, STAT4 and STAT6 were the most frequently mutated genes (Fig. 3A). There was a low to moderate correlation among members of the STAT family (Fig. 3B). Members of the STAT family play significant roles in cancer related pathways (Fig. 3C). In GBM, STATs were demonstrated to function in activating apoptosis pathway, EMT pathway, hormone ER pathway, RAS/MAPK pathway and RTK pathway, while inhibiting cell cycle pathway, DNA damage response pathway, hormone AR pathway and P13K/AKT pathway (Fig. 3C). GBM patients with high STAT5A expression leads to the sensitivity of 59 drugs or small molecules. Furthermore, GBM patients with high STAT5A expression leads to the sensitivity of 42 drugs or small molecules (Fig. 4).

OPEN IN VIEWER

Expression of STAT5A was high in GBM, making it a potential prognostic biomarker. Moreover, STAT5A was involved in activation and inhibition of key cancer related pathways and drug resistance. This implies that STAT5A plays significant roles in GBM pathogenesis and progression. Under this backdrop, STAT5A was selected for further analysis. Correlation between STAT5A and immune infiltrates in GBM was analyzed owing to the important roles the STAT family play in immune response [25, 26]. The levels of STAT5A were significantly correlated with abundance of B cell ( p = 1.91e-02), CD4 + T cells ( p = 9.55e-03), Macrophage ( p = 1.64e-04), Neutrophils ( p = 6.63e-04) and Dendritic cells ( p = 1.60e-16) (Fig. 5).

Expression of STAT5A was also strongly correlated with expression of gene biomarkers in GBM presented in Table 3. In reference to gene biomarkers of CD8+ T cell, data obtained from TIMER showed CD8B expression to be positively associated with STAT5A expression while data from CGGA showed the expression of both CD8A and CD8B to be positively associated with expression of STAT5A. Data obtained from both TIMER and CGGA indicated a strong correlation between expression of STAT5A and the expression of T cell (CD3D, CD3E, CD2), Monocyte (CD86 and CD115), TAM (CCL2, CD68, IL10) and M2 Macrophage (CD163, VSIG4, MS4A4A) gene biomarkers. From the analysis of gene biomarkers of Natural killer cell in TIMER, the levels of KIR2DL3, KIR2DL4, and KIR3DL1 were positively correlated with STAT5A expression. In both TIMER and CGGA, eight gene biomarkers of Dendritic cell were strongly associated with STAT5A expression. Similarly, the levels of STAT4, STAT1, TNF, GATA3, STAT6 and IL13 were positively correlatively with STAT5A expression. Two gene biomarkers (FOXP3 and TGFB1) of Treg and three gene biomarkers (PD-1, CTLA4 and GZMB) of exhibited a strong correlation with STAT5A expression. Therefore, STAT5A presents a suitable target for immunotherapy in GBM.

Figure 6A showed the genes significantly correlated with STAT5A in GBM. Among these genes, the top 50 most significant genes positively (Fig. 6B) and negatively (Fig. 6C) associated with STAT5A were also extracted. Moreover, Supplementary Fig. 1 showed the top three genes positively associated with STAT5A, which were TMEM199 (cor = 0.6939, p = 2.721e-23), FES (cor = 06883, p = 8.324e-23), and SLC24A6 (cor = 0.6829, p = 2.414e-22). Enrichment analyses of STAT5A, including GO function and KEGG pathways, were performed. The data of GO function analysis in Fig. 6D–F suggested that the functions of STAT5A in GBM mainly enriched in immune response-regulating signaling pathway, neutrophil and lymphocyte mediated immunity, regulation of immune effector process, immunological synapse, cytokine receptor binding, single-stranded DNA binding, and immunoglobulin binding. Furthermore, the result of KEGG pathway analysis suggested that the functions of STAT5A in GBM mainly enriched in cytokine-cytokine receptor interaction, NOD-like receptor signaling pathway, NF-kappa B signaling pathway, TNF signaling pathway, and Th17 cell differentiation (Fig. 6G and H).

OPEN IN VIEWER

Kinase and transcription factor target networks of STAT5A in GBM were also explored (Table 4). Here, kinases LCK, SYK, LYN, HCK, and TTK were identified as the five most significant kinase network targets. The PPI network showed kinase LCK target networks to be involved in regulation of immune response, antigen receptor-mediated signaling pathway, T cell activation and receptor signaling pathway (Fig. 7A). The five most significant transcription factor targets identified were V$ELF1_Q6, V$IRF_Q6, V$PEA3_Q6, V$PU1_Q6, and V$NERF_Q2. PPI network indicated that transcription factor ELF1 networks are involved in regulation of immune response, lymphocyte, leukocyte and B cell activation (Fig. 7B).

OPEN IN VIEWER